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INTRODUCTION 

The  mechanical  behavior  of  the  musculoskeletal  system  and 
the  associated  mathematical  models  have  become  an  important 
area  of  biomedical  research.  In  orthopaedics,  developments 
in  artificial  replacements  and  reconstructive  surgery  have 
demanded  greater  knowledge  of  the  mechanical  functions  of  the 
major  joints.  In  other  disciplines  such  as  sports  medicine, 
crash  protection,  and  Air  Force  related  applications  an 
increasing  interest  in  the  motions  and  forces  in  the  human 
body  is  noted.  This  need  for  a  better  understanding  of  the 
complicated  mechanical  behavior  of  biological  structures  has 
led  to  the  introduction  of  research  methods  and  mathematical 
tools  from  the  fields  of  life  sciences  as  well  as  applied  and 
theoretical  mechanics. 

In  attempting  to  understand  the  biodynamic  response  of 
the  human  body  subjected  to  expected  and/or  unexpected 
external  load  conditions,  properly  developed  mathematical 
models  can  provide  a  sound  basis  for  the  design  of  support- 
restraint  systems  and  vehicles  as  well.  The  most  sophisticat¬ 
ed  versions  of  these  mathematical  models  are  the  articulated 
and  multisegmented  total-human-body  models  which  initially 
appeared  in  the  crash  victim  simulation  literature.  These 
models  simulate  all  the  major  articulating  joints  and  segments 
of  the  human  body.  Representative  references  are  McHenry 
[1963],  Bartz  and  Butler  [1972],  Huston,  Hessel  and  Passerelo 
[1974],  Fleck,  Butler  and  Vogel  [1975]  and  Fleck  [1975], 

Currently,  the  Aerospace  Medical  Research  Laboratory 
(AMRL)  is  in  possession  of  an  Articulated  Total  Body  (ATB) 
model  which  is  maintained  and  used  in  Air  Force  related  appli¬ 
cations,  in  particular,  to  study  the  pilot  ejection  problem. 
The  ATB  model  has  some  special  features  which  allow  the 
capability  of  prescribing  a)  time-dependent  forces  on  the  body 
segments  to  simulate  aerodynamic  forces,  b)  joint  torques 
which  are  functions  of  both  flexure  and  azimuth  angles,  and 


c)  a  generalized  restraint  belt  subroutine  which  describes 
and  simulates  a  typical  Air  Force  harness  system. 

The  effectiveness  of  these  multisegmented  models  to 
accurately  predict  in-vivo  response  depends  upon  the  bio¬ 
mechanical  description  and  simulation  of  the  articulating 
joints.  This  study,  therefore,  concerns  with  the  analysis  of 
the  mechanical  behavior  of  the  major  articulating  joints  and 
the  development  of  mathematical  models  simulating  their 
dynamic  behavior.  In  general  mathematical  models  are  based  on 
physical  principles  and  consist  of  a  set  of  mathematical  rela¬ 
tions  among  relevant  parameters  of  the  system.  Assumptions 
and  simplifications  are  introduced,  ignoring  elements  assumed 
not  to  be  relevant  to  the  system's  behavior.  The  relevant 
parameters  of  the  system  are  defined  and  the  relations  are 
specified.  This  descriptive  model  is  then  expressed  as  a  set 
of  mathematical  relations  among  the  chosen  parameters.  The 
values  assigned  to  the  system's  parameters  are  selected  from 
the  literature  or  determined  by  measurement.  In  some  cases, 
the  parameters  cannot  be  measured  and  their  values  must  be 
estimated. 

Validation  of  a  model  is  established  when  the  model 
predictions  correlate  acceptably  with  data  in  the  literature 
and,  if  available,  with  results  of  experiments.  Within  the 
framework  of  the  present  study,  no  validation  experiments  were 
conducted,  so  that  model  predictions  could  be  compared  only 
with  experiments  reported  in  the  literature.  Considering  the 
conditions  of  the  reported  experiments  are  only  partly  known 
and  owing  to  the  variability  between  specimens,  such  a 
comparison  should  be  viewed  as  an  approximate  one. 

In  this  report,  a  rather  extended  discussion  of  the 
articulations  and  anatomical  descriptions  of  the  elbow, 
shoulder,  hip,  knee  and  ankle  joints  will  be  first  presented, 
with  special  emphasis  on  the  location  and  functional  aspects 
of  the  major  ligaments  of  each  joint.  This  is  followed  by  a 


description  of  the  articulating  surfaces  and  the  development 
of  a  measurement  technique  for  the  determination  of  articulat¬ 
ing  surface  equations  for  the  elbow,  hip,  knee  and  ankle 
joints.  Next,  a  constitutive  equation  representing  ligament 
characteristics  and  behavior  is  presented  and  the  attachment 
sites  of  the  ligaments  of  the  elbow,  hip,  knee  and  ankle  joints 
are  provided. 

General  two-  and  three-dimensional  mathematical  dynamic 
models  of  an  articulating  joint  are  then  developed  to  deter¬ 
mine  the  nature  of  motions  and  forces  between  two  body 
segments.  The  governing  equations  for  these  models  are  set 
of  highly  nonlinear  equations  and  their  numerical  solutions 
are  discussed  in  some  detail.  This  is  followed  by  a  specific 
application  to  a  two-dimensional  dynamic  model  of  the  human 
knee  joint.  The  numerical  results  from  this  model  are 
presented  to  illustrate  the  effects  of  duration  and  shape  of 
the  dynamically  applied  loads  on  the  response  of  the  joint. 
Special  attention  has  been  given  to  the  ligament  and  contact 
forces,  the  location  of  contact  points,  anterior-posterior 
displacements  and  the  comparison  between  the  internal  and  the 
external  energy  of  the  system.  The  results  are  compared  with 
experimental  data  from  the  literature  and  the  validation  of 
the  model  is  established.  The  report  is  concluded  with  a 
discussion  of  extensions  of  the  model  and  its  possible  impli¬ 
cations  on  future  research. 

ARTICULATION  AND  ANATOMICAL  DESCRIPTION  OF 

ELBOW,  SHOULDER,  HIP,  KNEE  AND  ANKLE  JOINTS 

Realistic,  accurate  mathematical  modelling  of  the  major 
articulating  joints  of  the  human  body  requires  a  comprehensive 
knowledge  and  understanding  of  the  physical  behavior  and  the 
anatomical  characteristics  of  each  joint.  In  a  previous 
report  [Engin,  1979a]  a  survey  was  provided  for  various  major 
human  joint  models  including  a  single  degree  of  freedom  hinge 
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or  revolute  joint,  a  spherical  joint  limited  to  two  degrees  of 
freedom,  three  degrees  of  freedom  planer  joint,  three  degrees 
of  freedom  ball  and  socket  joint,  and  a  general  six  degrees  of 
freedom.  Passive  and  active  force  and  moment  response  of 
major  human  joints,  associated  torques  about  the  long-bone 
axes  of  these  joints  and  some  aspects  of  joint  modeling  were 
reported  in  the  literature  by  the  senior  author  in  a  series  of 
articles  [Engin,  1979b;  Engin  et  al .  1979c;  Engin,  1979d; 

Engin  and  Kaleps,  1980a;  Engin,  1980b;  Engin  and  Peindl,  1980c; 
Engin,  Akkas  and  Kaleps,  1980d;  Engin,  1981a§b;  Engin  and 
Moeinzadeh,  1981c].  The  research  works  presented  in  these 
articles  were  performed  with  some  obvious  limitations  on  live 
subjects  by  means  of  specially  designed  experimental  apparatus. 

In  the  following  paragraphs,  descriptions  of  the  essential 
anatomical  and  functional  aspects  of  the  elbow,  shoulder,  hip, 
knee  and  ankle  joints  will  be  presented.  For  each  joint,  the 
physical  structure  and  the  movements  of  the  articulating  seg¬ 
ments  are  described  and  the  ligaments  having  a  significant 
contribution  to  the  integrity  and  function  of  each  joint  are 
defined.  A  number  of  illustrative  figures  are  presented  for 
each  joint  with  segments  appropriately  identified  by  the 
commonly  accepted  medical  terminology.  The  anatomical  and 
functional  descriptions  provided  below  were  taken  from  various 
sources  such  as  Gray  [1973],  Grant  [1962],  and  Wells  [1971]. 

ELBOW  JOINT 

The  elbow  joint  is  a  uni-axial  (hinge)  joint.  It  is 
composed  of,  proximally,  the  trochlea  and  capitulum  of  the 
humerus,  and  distally,  the  trochlear  notch  of  the  ulna  and  the 
head  of  the  radius.  The  trochlea  of  the  humerus  is  convex, 
anteroposterior ly ,  and  concave,  side  to  side,  fitting  into  the 
trochlear  notch  of  the  ulna.  The  spherical  capitulum  of  the 
humerus  fits  into  the  concave  head  of  the  radius.  The  cavity 
of  the  elbow  joint  is  continuous  with  the  superior  radio-ulnar 
joint  where  the  head  of  the  radius  fits  into  the  radial  notch 
of  the  ulna. 


Movements  of  the  elbow  consist  of  flexion  and  extension, 
determined  by  the  shape  of  the  trochlear  surfaces:  the  medial 
portion  of  the  trochlea  projects,  distally,  further  than  the 
lateral  portion.  Flexion  is  limited  by  the  soft  tissues  of 
the  arm,  whereas  extension  is  limited  by  the  olecranon  of  the 
ulna  contacting  the  base  of  the  olecranon  fossa  (Figure  1) . 

The  capsule  of  the  elbow  joint  is  thin  and  covered  by 
muscles,  attaching  anteriorly  to  the  humerus  slightly  above 
the  coronoia  process  and  to  the  annular  ligament  around  the 
head  of  the  radius.  Posteriorly,  attachment  is  slightly  above 
the  capitulum  of  the  humerus,  to  the  olecranon  fossa,  the 
olecranon  upper  margins  and  to  the  capsule  of  the  superior 
radio-ulnar  joint. 


The  fibrous  capsule  is  lined  by  a  synovial  membrane  which 
extends  into  the  coronoid,  olecranon  and  radial  fossae  and 
between  the  ulna  and  radius,  being  continuous  with  the  syno¬ 
vial  membrane  of  the  superior  radio-ulnar  joint.  Pads  of  fat 
are  present  between  the  fibrous  capsule  and  the  synovial 
membrane.  These  fat  pads  fill  the  fossae:  the  olecranon  fossa 
during  flexion,  and  the  coronoid  and  radial  fossa  during 
extension . 

The  capsule  is  strengthened  by  three  ligaments: 

1.  ulnar  collateral  ligament  (medial  ligament); 

2.  radial  collateral  ligament  (lateral  ligament);  and 

3.  annular  ligament. 

Ulnar  Collateral  Ligament 

The  ulnar  collateral  ligament  is  a  roughly  triangular 
thick  band,  composed  of  a  strong  anterior  band  and  a  weaker 
middle  and  transverse  sheet  (Figure  2) .  It  extends  from  the 
medial  epicondyle  of  the  humerus  to  an  attachment  along  the 
coronoid  process  and  the  olecranon  of  the  ulna.  The  anterior 
portion  is  nearly  a  cord,  being  taut  in  extension.  The  poste¬ 
rior  portion  attaches  to  the  distal  and  posterior  of  the 
medial  epicondyle  and  to  the  medial  margin  of  the  olecranon. 
This  portion  of  the  ligament  is  a  weaker  sheet  which  is  taut 
in  flexion.  An  oblique  band  extends  between  the  olecranon 
and  the  coronoid  process,  deepening  the  socket  for  the 
trochlea  of  the  humerus. 


Radial  Collateral  Ligament 

This  ligament  is  attached  to  the  lateral  epicondyle  of 
the  humerus,  to  the  trochlear  notch  of  the  ulna  and  to  the 
annular  ligament.  It  does  not  attach  directly  to  the  radius 
so  that  rotation  of  the  radius  is  permitted  in  pronation  and 
supination  of  the  forearm  (Figure  3). 


Figure  3.  Lateral  aspect  of  the  right  elbow  joint 


Annular  Ligament 

This  strong  band  encircles  the  head  and  neck  of  the 
radius,  attaching  to  the  anterior  and  posterior  margins  of 
the  radial  notch  of  the  ulna.  The  radial  collateral  ligament 
blends  into  the  proximal  margin  of  the  annular  ligament 
(Figure  3).  With  the  oblique  cord,  which  extends  medially  and 
upwards  from  the  radial  tuberosity  to  the  coronoid  process  of 
the  ulna  (Figure  3),  the  annular  ligament  maintains  the  radius 
head  close  to  the  radial  notch  of  the  ulna. 


SHOULDER  COMPLEX 

The  shoulder  complex  is  composed  of  four  independent 
articulations  among  the  bones  of  the  complex:  the  clavicle, 
scapula,  humerus  and  the  thorax  (Figure  4).  The  shoulder 
girdle  is  composed  of  the  clavicle  and  scapula.  There  are  two 
clavicular  articulations:  the  sternoclavicular  joint,  where 
the  clavicle  articulates  with  the  manubrium  of  the  sternum,  and 
the  acromioclavicular  joint,  where  the  clavicle  articulates 
with  the  acromion  process  of  the  scapula.  The  glenohumeral 
joint  is  a  ball  and  socket  joint  composed  of  the  humerus  and 
the  glenoid  cavity  of  the  scapula.  The  final  articulation  is 
not,  per  se,  a  joint  but  is  the  scapulothoracic  articulation 
of  the  scapula  over  the  thorax. 

Sternoclavicular  Joint 

This  is  a  saddle-type  joint  with  both  concave  and  convex 
curvatures.  The  proximal  end  of  the  clavicle  is  separated 
from  the  manubrium  of  the  sternum  by  a  constant  thickness, 
intra-articular  meniscus.  The  fibrous  joint  capsule  is 
strengthened  by  the  anterior  and  posterior  sternoclavicular 
ligaments.  Further,  both  left  and  right  clavicles  are  joined 
by  the  interclavicular  ligament  running  over  the  sternal 
notch  (Figure  4) .  The  inferior  portion  of  the  clavicle  con¬ 
nects  to  the  first  costal  cartilage  at  the  costal  tuberosity 
by  means  of  the  costoclavicular  ligament.  The  sternoclavicular 
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Figure  4.  Anterior  view  of  the  right  shoulder  complex. 


Figure  5.  Posterior  aspect  of  the  right  shoulder  complex 


joint  possesses  three  degrees  of  freedom  with  axes  in  the 
sagittal  and  frontal  planes  and  the  bone-axes  of  the  clavicle. 
Elevation  of  the  clavicle  is  limited  by  the  lower  portion  of 
the  joint  capsule  and  the  costoclavicular  ligament.  Depression 
is  limited  by  the  upper  portion  of  the  joint  capsule  and  the 
interclavicular  ligament. 

Acromioclavicular  Joint 

This  articulation  between  the  distal  end  of  the  clavicle 
and  the  acromion  of  the  scapula  is  surrounded  by  a  fibrous 
capsule.  The  articulation  provides  little  stability  so  that 
the  ligamentous  structures  are  the  primary  stabilizers. 
Reinforcing  the  capsule  are  the  superior  and  inferior  acromio¬ 
clavicular  ligaments  (Figure  4).  Additionally,  the  clavicle 
connects  to  the  scapula  by  the  conoid  and  trapezoid  portions 
of  the  coracoclavicular  ligament  and  by  the  coracoacromial 
ligament  (Figure  5) .  The  proximity  of  the  coracoid  process 
of  the  scapula  to  the  clavicle  (and  the  possible  cartilagenous 
formation  between  them) ,  sometimes  is  referred  to  as  the 
coracoclavicular  joint  with  the  entire  region  being  called  the 
claviscapular  joint. 

Glenohumeral  Joint 

This  ball  and  socket  joint  between  the  head  of  the 
humerus  and  the  glenoid  fossa  of  the  scapula  is  surrounded  by 
a  loose  sleeve  composed  of  the  joint  capsule  and  its  capsular 
ligaments.  The  contact  surface  is  remarkably  small  with  the 
head  of  the  humerus  possessing  a  much  greater  articulating 
surface  than  the  glenoid  fossa  of  the  scapula.  In  addition  to 
the  gelnohumeral  capsular  ligaments,  the  coracohumeral  liga¬ 
ment  over  the  superior  aspect  of  the  joint  provides  reinforce¬ 
ment  as  the  humerus  is  suspended  along  side  of  the  torso  as 
well  as  checks  outward  rotation  of  the  humerus  (Figure  6) . 
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Figure  6.  Superior  aspect  of  the  right  shoulder  complex. 


HIP  JOINT 

The  hip  joint  is  a  multi-axial,  ball  and  socket  joint, 
and  therefore  possesses  three  degrees  of  freedom.  For  conve¬ 
nience,  movements  can  be  considered  to  be  about  three  mutually 
perpendicular  axes:  transverse,  anteroposterior  and 
longitudinal.  Flexion  and  extension  occur  about  the  trans¬ 
verse  axis;  flexion  being  forward  movement  and  extension  being 
backward.  Flexion  is  limited  by  tension  on  the  hamstrings  and 
soft  tissues;  extension  is  limited  by  the  iliofemoral  and 
pubofemoral  ligaments.  Abduction,  movement  of  the  thigh  away 
from  the  midline  of  the  body,  and  adduction,  movement  toward 
the  midline,  occurs  about  the  anteroposterior  axis.  Abduction 
is  limited  by  tension  of  the  adductor  muscles  and  by  contact 
of  the  greater  trochanter  with  the  acetabulum  (Figure  7) . 
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Adduction  is  limited  by  the  opposing  leg  and  by  the  iliofemoral 
and  ischiofemoral  ligaments.  Lateral  and  medial  rotation, 
movement  of  the  anterior  surface  of  the  thigh  laterally  and 
medially,  respectively,  occurs  about  the  longitudinal  axis. 
Lateral  rotation  is  limited  by  the  iliofemoral  and  pubofemoral 
ligaments  (Figure  8);  medial  rotation,  by  the  ischiofemoral 
and  iliofemoral  ligaments  (Figure  10)  . 

Analysis  of  these  movements  must  take  into  account  the 
length  and  angulation  of  the  neck  of  the  femur  in  relation  to 
the  long  axis  of  the  femoral  shaft.  In  flexion  or  extension, 
the  femoral  head  rotates  about  a  transverse  axis  within  the 
acetabulum.  Medial  and  lateral  rotations  occur  about  a  longi¬ 
tudinal  axis  through  the  head  of  the  femur  and  the  lateral 
condyle,  when  the  foot  is  weight-bearing.  Thus,  the  medial 
condyle  moves  posteriorly  and  the  greater  trochanter  moves 
anteriorly  in  relation  to  this  axis  during  medial  rotation 
However,  when  the  foot  is  free  or  not  weight-bearing,  rotation 
may  occur  about  variable  axes  through  the  federal  head. 

Abduction  and  adduction  are  produced  about  .i.i  anteroposterior 
axis  through  the  approximate  center  of  the  femoral  head. 

Because  of  shifting  axes,  it  is  sometimes  convenient  to 
view  movements  as  occurring  about  mechanical  axes  through  the 
femoral  neck  and  the  approximate  center  of  the  femoral  head. 
Thus,  extension  and  flexion  can  be  viewed  as  spins,  producing 
respectively,  tauting  (spiralizing)  and  relaxing  (straightening) 
of  the  ligaments  and  the  capsule. 

The  joint  itself  is  composed  of  the  head  of  the  femur  and 
the  acetabulum  of  the  hip.  The  femoral  head  is  approximately 
two-thirds  of  a  sphere.  The  acetabulum  is  shaped  like  a 
horseshoe,  closed  by  a  non-articulating  pad  of  fat.  The 
socket  is  deepened  by  the  acetabular  Inbrum,  a  fibrocartilag¬ 
inous  rim,  and  the  transverse  ligament,  which  completes  the 
encapsulation  of  the  head  of  the  femur.  The  femoral  head 
aligns  obliquely  upwards,  medially  and  slightly  forwards. 
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The  capsule  attaches  slightly  beyond  the  acetabular 
labrum,  blending  with  the  labrum  and  the  transverse  ligament 
anteriorly  and  inferiorly.  Femoral  attachments  are  along  the 
intertrochanteric  line,  anteriorly,  and  to  the  neck,  medially 
to  the  obturator  externus.  Three  ligaments  strengthen  the 
capsule : 

1.  iliofemoral  ligament; 

2.  pubofemoral  ligament;  and 

3.  ischiofemoral  ligament. 

Iliofemoral  Ligament 

Triangular  in  shape  and  of  great  strength,  the  iliofemoral 
ligament  covers  the  anterior  portion  of  the  joint.  It  attaches 
to  the  ilium  above  the  acetabular  rim  and  to  the  lower  portion 
of  the  anterior  inferior  iliac  spine  (Figure  9).  The  ligament 
broadens  and  diverges,  inferiorly,  to  form  two  bands.  The 
superior  or  lateral  band  attaches  to  the  upper  part  of  the 
intertrochanteric  line  and  is  sometimes  referred  to  as  the 
il iotrochanteric  ligament.  The  inferior  or  medial  band 
attaches  to  the  lower  portion  of  the  intertrochanteric  line. 

The  iliofemoral  ligament  is  referred  to  as  the  "Y"  ligament 
due  to  its  inverted  "Y”-shape.  This  ligament  checks  extension 
of  the  joint,  as  well  as  both  lateral  and  medial  rotation. 

Pubofemoral  Ligament 

The  pubofemoral  ligament  attaches  medially  to  the 
anterior  acetabular  rim  and  the  superior  pubic  ramus,  and 
crosses  to  the  inferior  of  the  neck  of  the  femur,  attaching  to 
the  top  of  the  lesser  trochanter  (Figure  8) .  Somewhat  trian¬ 
gular  in  shape,  this  ligament  checks  abduction,  extension  and 
medial  rotation. 

Ischiofemoral  Ligament 

This  ligament  is  less  differentiated  than  the  preceding. 

It  attaches  to  the  ischium  posteriorly  and  inferiorly  to  the 
acetabulum,  passing  over  the  superior  and  posterior  of  the 
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neck  of  the  femur  (Figure  10).  The  superior  portion  is 
generally  horizontal  and  the  posterior  or  lower  portion 
spirals.  The  ischiofemoral  ligament  limits  medial  rotation 
and  adduction. 

Teres  Femoris  (Head  of  Femur)  Ligament 

Structurally,  this  ligament  is  of  little  importance. 
Passing  from  the  deep  aspect  of  the  transverse  ligament  to  a 
depression  in  the  head  of  the  femur;  it  functions  primarily 
as  a  conduit  for  blood  vessels  (Figure  7) .  In  extreme 
abduction,  this  ligament  becomes  taut,  but  not  before  the 
iliofemoral  ligament  has  become  taut  and  possibly  not  until 
the  iliofemoral  ligament  fails. 

KNEE  JOINT 

A  condyloid,  synovial  joint,  the  knee  is  the  articulation 
of  the  distal  condylar  surfaces  of  the  femur,  the  proximal 
condylar  surfaces  of  the  tibia  and  the  posterior  surface  of 
the  patella.  While  the  primary  movement  is  hinge-like,  some 
rotation  does  occur.  Flexion,  backward  movement  of  the  thigh 
or  leg,  and  extension,  the  opposite  movement,  occur  about  a 
moving  transverse  axis.  This  axis  moves  backward  during 
flexion  due  to  the  curvatures  of  the  femoral  condyles. 

At  complete  flexion,  the  posterior  femoral  condylar 
surfaces  articulate  with  the  posterior  tibial  condylar  sur¬ 
faces  and  with  the  posterior  portions  of  the  menisci.  During 
extension  with  the  tibia  fixed,  the  femoral  condyles  roll 
forward  while  simultaneously  sliding  backward  on  the  tibial 
condylar  surfaces.  The  contact  areas  between  these  surfaces 
increases  as  the  curvature  of  the  femoral  condyles  decreases. 
Movement  on  the  lateral  condyle  ends  before  extension  is 
complete,  while  movement  on  the  medial  condyle  continues  since 
the  lateral  articular  surface  of  the  lateral  condyle  is  short¬ 
er  than  the  medial.  This  continued  movement  of  the  medial 
condyle  causes  the  femur  to  rotate  medially  about  a 
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longitudinal  axis  through  the  lateral  condyle,  causing  the 
collateral  and  popliteus  oblique  ligaments  to  become  taut. 
Thus,  lateral  rotation  becomes  a  precursor  to  flexion. 

During  extension,  there  is  a  lateral  rotation  of  the 
tibia  with  respect  to  the  femur.  Conversely,  medial  rotation 
of  the  tibia  occurs  in  the  beginning  of  flexion.  When 
approaching  complete  extension,  the  anterior  portions  of  the 
menisci  are  pushed  forward  by  the  femur  and  become  less 
curved.  The  opposite  occurs  in  flexion.  As  the  tibial  col¬ 
lateral  ligament  becomes  taut  during  extension,  it  pulls  the 
medial  meniscus  outward.  The  lateral  meniscus  is  drawn  out¬ 
ward,  away  from  the  femoral  and  tibial  condyles,  by  the 
popliteus  which  is  attached  to  the  posterior  of  the  lateral 
meniscus . 

The  cruciate  ligaments  are  taut  in  most  positions  of  the 
knee  and  prevent  anteroposterior  displacement  of  the  tibia  in 
relation  to  the  femur.  During  rotational  movements,  the 
cruciates  twist  and  untwist  around  each  other.  In  full 
flexion,  the  anterior  cruciate  ligament  is  relaxed  and  in  full 
extension  the  posterior  cruciate  is  relaxed.  The  collateral 
ligaments  are  relaxed  when  the  knee  is  flexed  to  a  ninety- 
degree  angle,  thus  allowing  rotation  about  a  vertical  axis. 

The  articular  surfaces  of  the  femoral  condyles  are  convex 
anteroposterior ly  and  from  s ide - to- s ide  ,  being  more  marked  in 
the  posterior  portion  of  the  anteroposterior  curvature.  The 
tibial  surfaces  arc  comparatively  flat,  being  deepened  by  the 
wedge-shaped  menisci.  The  patellar  surface  of  the  femur  is 
also  convex  from  s ide - to - s ide  ,  with  the  lateral  condyle 
extending  further  forward  and  upward. 

The  medial  and  lateral  menisci  are  two  crescent  -  shaped 
fibrocartilaginous  structures  attached  to  the  upper  surface 
of  the  tibia  by  the  coronary  ligaments.  The  inferior  surface 
is  flattened  while  the  superior  surface  is  concave,  thus 
deepening  the  sockets  for  the  femoral  condyles.  The  outside 
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edges  are  firmly  attached  to  the  tibia  through  the  capsule  to 
the  tibial  condyles.  The  inner  edges  are  thin  and  free.  The 
anterior  ends  are  attached  by  the  transverse  ligament.  The 
posterior  portion  of  the  lateral  meniscus  (Figure  11)  is  con¬ 
tinuous  with  the  posterior  cruciate  ligament  and  attaches  to 
the  popliteous  tendon.  The  medial  meniscus  attaches  to  the 
tibial  collateral  ligament  (Figure  14) .  The  lateral  is  broader 
and  its  ends  closer  than  the  medial  meniscus.  The  lateral  is 
more  nearly  circular  while  the  medial  meniscus  is  more 
elliptical . 

Four  major  ligaments  lend  stability  to  the  knee  joint. 

1.  tibial  collateral  (medial)  ligament; 

2.  fibular  collateral  (lateral)  ligament; 

3.  anterior  cruciate  ligament;  and 

4.  posterior  cruciate  ligament. 

Tibial  Collateral  (Medial)  Ligament 

This  is  a  broad,  flat  band  in  the  medial  portion  of  the 
capsule.  It  is  attached  proximately  to  the  medial  epicondyle 
of  the  femur  below  the  adductor  tubercle,  and  it  broadens  to 
an  attachment  at  the  medial  condyle  and  upper  body  of  the 
tibia  (Figure  11).  Its  deep  fibers  attach  to  the  medial 
meniscus  periphery. 


Fibular  Collateral  (Lateral)  Ligament 

Unlike  the  tibial  collateral  ligament,  the  fibular 
collateral  ligament  is  distinctly  separate  from  the  fibrous 
capsule  (Figure  11).  It  is  a  strong,  rounded  cord,  attached 
to  the  lateral  epicondyle  above  the  groove  for  the  popliteous 
tendon  and  passes  to  the  lateral  side  of  the  head  of  the 
fibula.  The  popliteous  tendon  lies  below  it,  separating  it 
from  the  lateral  meniscus. 

Anterior  Cruciate  Ligament 

The  anterior  cruciate  ligament  attaches  to  the  tibia 
anterior  to  the  intercondylar  eminence,  between  the  menisci. 
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partly  blending  with  the  anterior  end  of  the  lateral  meniscus 
(Figure  12).  It  crosses  upward,  backward  and  laterally, 
twisting  on  itself,  fanning  out  to  attach  to  the  medial  aspect 
of  the  femoral  lateral  condyle.  It  is  anterolateral  to  the 
posterior  cruciate  ligament.  It  holds  the  femur  from  sliding 
backward,  prevents  hyperextension  of  the  knee  and  checks 
medial  rotation  of  the  femur  when  the  leg  is  weight-bearing. 

Posterior  Cruciate  Ligament 

The  posterior  cruciate  ligament  attaches  to  the  posterior 
intercondylar  area  of  the  tibia  and  the  lateral  meniscus, 
crossing  upward,  forward  and  medially  to  attach  to  the  lateral 
portion  of  the  femoral  medial  condyle  (Figure  13).  It  is 
stronger,  shorter  and  less  oblique  in  direction  than  the 


Figure  13.  Posterior  aspect  of  the  right  knee  joint. 


anterior  cruciate  ligament.  It  holds  the  femur  from  sliding 
forward.  The  two  cruciates  are  slightly  twisted  around  each 
other  except  when  the  knee  is  fully  extended.  At  this 
position,  the  tibia  is  laterally  rotated  in  relation  to  the 
femur . 

ANKLE  JOINT 

The  ankle  or  talocrural  joint  is  a  uni-axial  hinge  joint 
formed  by  the  articulation  of  the  talus  with  a  three-sided 
socket  composed  of  the  distal  surface  of  the  tibia  and  the 
articular  surfaces  of  the  tibial  and  fibular  malleoli  with  the 
inferior  transverse  tibiofibular  ligament,  posteriorly.  The 
tibia  and  fibula  are  firmly  united  at  the  inferior  tibio¬ 
fibular  joint. 

The  hinge- like  movement  of  the  ankle  occur  about  an  axis 
through  the  body  of  the  talus  (Figure  14),  which  is  slightly 
oblique,  passing  forward,  medial  to  lateral.  Dorsiflexion 
(extension)  is  raising  the  forepart  of  the  foot  while  plantar 
flexion  (flexion)  is  the  lowering  of  the  foot.  There  is 
maximal  congruence  of  the  joint  surfaces  and  maximal 
ligamentous  tension  in  dorsiflexion. 

The  joint  capsule  is  thin  anteriorly  and  posteriorly  with 
lateral  ligaments.  There  are  deep  fatty  pads  in  the  anterior 
and  posterior  portions  of  the  joint.  The  joint  cavity  extends 
upward  between  the  tibia  fibula  for  a  few  millimeters.  The 
anterior  portion  of  the  capsule  is  attached  to  the  tibia  near 
the  articular  surface  and  to  the  neck  of  the  talus  near  its 
head . 

The  integrity  of  the  ankle  is  determined  in  part  by  the 
bony  structure  and  in  part  by  two  ligaments: 

1.  medial  ligament  and 

2.  lateral  ligament. 

Medial  Ligament 

The  medial  ligament  (deltoid  ligament)  is  a  strong,  thick 
triangular  band  connecting  the  medial  malleolus  of  the  tibia 
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to  several  tarsal  bones  (Figure  14) .  The  deep  portion  passes 
from  the  malleolus  to  the  medial  surface  of  the  talus.  It  is 
also  known  as  the  anterior  tibiotalar  ligament.  The  superfi¬ 
cial  portion  consists  of  three  parts:  (1)  the  anterior  or 
tibionavicular  portion  passes  to  the  tuberosity  of  the  navic¬ 
ular  bone,  blending  with  the  medial  margin  of  the  plantar 
calcaneonavicular  ligament;  (2)  the  middle  or  tibiocalcanean 
portion  crosses  to  the  sustentaculum  cali  of  the  calcaneous; 
and  (3)  the  posterior  or  posterior  tibiotalar  portion  attaches 
to  the  medial  tubercle  of  the  talus. 

Lateral  Ligament 

The  lateral  ligament  consists  of  three  parts.  (1)  The 
anterior  talofibular  ligament  attaches  to  the  anterior  margin 
of  the  lateral  malleolus  of  the  fibula,  and  to  the  anterior 
of  the  talus  between  its  neck  and  articular  surface  (fibular 
articulation)  (Figure  15) .  (2)  The  calcaneof ibular  ligament 

passes  down  and  back  from  the  tip  of  the  lateral  malleolus  to 
the  lateral  aspect  of  the  calcaneous,  posterior  to  the 
peroneal  tubercle.  (3)  The  posterior  talofibular  ligament 
passes  from  the  malleolar  fossa  to  the  lateral  tubercle  of  the 
talus. 

GEOMETRY  OF  THE  ARTICULATING  SURFACES 

Articulating  surfaces  play  a  major  role  in  the  physical 
motion  of  a  joint.  It  is  assumed  that  the  deformations  in 
the  articular  surfaces  do  not  affect  relative  motions  and 
forces  in  the  joint.  This  assumption  is  based  on  the  consid¬ 
eration  that  the  deformations  of  the  cartilage  layer  caused  by 
the  contact  between  the  articulating  bones  is  relatively 
slight  compared  to  the  range  of  motions  in  the  joint  (Wismans 
[1980]).  Therefore,  the  articular  surfaces  are  represented 
by  rigid  surfaces  and  the  contact  areas  between  the  articulat¬ 
ing  surfaces  are  reduced  to  contact  points. 
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No  accurate  quantitative  data  for  the  geometry  of  the 
articulating  surfaces  of  the  human  joints  are  available  in 
the  literature.  The  knee  joint  is  relatively  the  most  studied 
joint  and  some  techniques  have  been  developed  for  measurements 
of  its  articular  surfaces.  Seedham,  et  al.,  [1972a]  studied 
the  femoral  and  tibial  condyles  by  making  a  plastic  mold  of 
the  condyles.  These  mouldings  were  cut  in  the  sagittal 
planes,  resulting  in  a  number  of  contours  of  the  articular 
surface.  In  the  work  of  Wismans,  et  al.,  [1980],  using  dial 
gages,  three-dimensional  coordinates  of  50-200  points  on  each 
condyle  of  the  knee  were  determined  and  the  surfaces  were 
approximated  by  mathematical  functions  representing  the  col¬ 
lected  data  points.  Most  other  studies  in  this  field  are 
restricted  to  the  determination  of  a  number  of  rough  dimen¬ 
sions  from  roetgen  photographs  (Erkman  and  Walker  [1974]  , 
See-dham,  et  al.  [1972b]).  Several  other  measuring  techniques 
are  also  given  in  a  survey  by  Wismans  and  Struben  [1977]  and 
in  Devens  [1979],  where  special  attention  is  paid  to  optical 
methods . 

In  this  study,  coordinates  of  a  large  number  of  points 
on  each  of  the  articular  surfaces  of  the  elbow,  hip,  knee  and 
ankle  joints  are  determined  using  a  sonic  digitizing 
technique.  This  technique,  its  measuring  apparatus  and  pro¬ 
cedures  will  be  presented  below.  A  procedure  for  approximat¬ 
ing  the  articular  surfaces  by  mathematical  functions  will  then 
be  discussed  and  a  brief  description  of  the  articulating 
surfaces,  the  location  and  orientation  of  coordinate  systems 
and  the  mathematical  functions  representing  the  elbow,  hip, 
knee  and  ankle  joint  surfaces  will  be  presented. 

MEASURING  TECHNIQUE 

Coordinates  of  a  large  number  of  points  on  the 
articulating  surfaces  are  determined  using  a  Graf/Pen  Sonic 
Digitizer.  Sonic  digitizing  is  the  process  of  converting 
information  on  location  or  position  in  one,  two  or  three 
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dimensions,  to  digital  values  suitable  for  data  processing, 
storage  or  transmission.  The  system  used  to  accomplish  this 
conversion  consists  of  a  stylus,  two  or  three  microphone/ 
sensor  assemblies  (for  two-  and  three-dimensional  data 
conversion,  respectively),  an  electronic  control  unit  and  a 
generator/multiplexer  unit  which  is  used  to  select  and  power 
sonic  impulse  emitters.  As  an  example  of  the  digitizer's 
operation,  let  us  first  consider  the  two-dimensional  mode  of 
operation. 

The  two-dimensional  sensor  assembly  consists  of  two 
perpendicular,  linear  microphones  as  shown  in  Figure  16. 

These  sensors  define  a  planar  effective  working  area  of  ap¬ 
proximately  35  cm  x  35  cm.  The  Graf/Pen  uses  impulses  gener¬ 
ated  at  the  tip  of  the  emitter  to  calculate  its  position  in 


Figure  16.  The  two-dimensional  microphone/sensor  assembly  of 
the  Graf/Pen  sonic  digitizer. 
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the  working  area.  The  times  required  for  the  sound  waves  to 
reach  the  two  microphone/ sensors  are  converted  into  distance 
measurements  (the  x  and  y  coordinates) .  The  measurements  are 
then  transmitted  to  the  computer  in  binary  coded  decimal 
(BCD)  Cartesian  form.  The  two-dimensional  mode  is  primarily 
used  to  digitize  joint  surface  curves  which  are  obtained  from 
x-ray  projections  (see  Figure  32)  or  physical  cross-sections 
of  the  joint.  Additionally,  the  two-dimensional  microphone/ 
sensor  assembly  provides  a  menu  capability  of  alphanumeric 
data  entry  to  the  computer. 

An  obvious  advantage  of  using  sound  as  a  ranging  device 
is  that  digitization  need  not  be  confined  to  a  plane.  The 
three-dimensional  microphone/sensor  unit  utilized  for  the 
digitization  of  the  articulating  surfaces  consists  of  four 
linear  microphones  arranged  in  a  planar,  rectangular  manner. 
These  microphones  define  a  three-dimensional  effective  work¬ 
ing  volume  of  approximately  150  cm  x  75  cm  x  180  cm,  along 
the  x,  y  and  z  directions,  respectively.  The  three-dimensional 
mode  of  operation  is  similar  to  that  of  the  two-dimensional 
mode.  In  the  three-dimensional  set-up,  however,  the  dis¬ 
tances  measured  are  slant  ranges  to  each  of  the  coplanar 
sensors.  Thus,  the  information  generated  by  each  sensor 
represents  the  radius  of  a  circular  arc  which  includes  the 
impulse  source  {the  tip  of  the  sonic  emitter)  and  is  in  a 
plane  perpendicular  to  the  sensor.  Four  sensors  are  used 
merely  for  the  purposes  of  accuracy.  The  digitizer  examines 
the  signals  from  all  four  sensors,  selects  the  three  smallest 
signals  and  disregards  the  fourth.  The  location  of  the  sonic 
emitter  is  then  calculated  as  being  at  the  intersection  of 
the  three  smallest  arcs.  The  three  slant  ranges  are  easily 
converted  into  Cartesian  x,  y  and  z  coordinates  by  a  micro¬ 
processor  in  the  control  unit  and  converted  for  transmission. 

Using  the  joints  of  a  full-size  human  skeleton,  each 
articulating  joint  segment  was  placed  on  the  apparatus  shown 


in  Figure  17,  with  its  articular  surface  facing  the 
microphone/sensor  unit.  Digitizing  each  surface,  the  coordi¬ 
nates  of  a  large  number  of  points  (50-200)  on  the  surfaces 
were  measured  with  respect  to  the  microphone/sensor  unit's 
coordinate  system.  The  output  was  recorded  on  an  LA120 
terminal,  manufactured  by  Digital  Equipment  Corporation. 
Applying  simple  mathematics,  these  coordinates  were  then 
transformed  into  the  local  coordinate  system  of  the  particular 
joint  segment  under  study. 


Figure  17.  Experimental  set-up  for  measuring  the  three 
dimensional  geometry  of  the  articulating  surfaces. 


MATHEMATICAL  DESCRIPTION  OF  THE  ARTICULATING  SURFACES 

Based  on  the  above  procedure,  coordinates  (Xq ,  yg ,  z^)  of 
a  suffici  it  number  of  points  on  each  of  the  articulating 
surfaces  are  determined  in  their  respective,  local  coordinate 
systems.  The  aim  is  to  find  a  realistic  and  simple  represen¬ 
tation  in  the  form  of  a  mathematical  approximation  for  the 
geometry  of  these  surfaces.  In  this  study,  a  polynomial  of 
degree  n  (n>l)  is  used  and  expressed  as: 

y  =  f(x,z)  =  l  f  a xpzq  (1) 

p=0  q=0  Pq 

The  coefficients  ap^  of  equation  (1)  have  to  be  determined  in 
such  a  way  that  for  each  point  Q,  the  coordinate  yg  is 
approximated  as  accurately  as  possible  by  y(Xg,Eg).  In  the 
present  work,  the  coefficients  ap^  are  obtained  by  means  of 
statistical  operations  using  the  79. 3A  version  of  the  GLM 
procedure  of  the  Statistical  Analysis  System  (SAS)  subroutines 
(SAS  [1979]).  The  GLM  procedure  uses  the  principle  of  least 
squares  (Goodnight  and  Harvey  [1978],  Draper  and  Smith  [1966], 
Graybill  [1961])  to  fit  linear  models  and  provides  an  output 
data  set  containing  (a)  predicted  and  residual  values  from  the 
analysis,  (b)  standard  deviation  and  (c)  percentage  accuracy 
of  the  fit.  Summarized  values  for  the  degree  of  the 
polynomial,  n,  the  standard  deviation,  o,  and  the  percentage 
accuracy,  R,  of  the  models  obtained  for  the  articulating  sur¬ 
faces  of  the  elbow,  hip,  knee  and  ankle  joints  are  presented 
in  Table  1. 

ARTICULATING  SURFACES  OF  THE  ELBOW  JOINT 
Anatomical  Description 

The  two  main  articulating  surfaces  of  the  elbow  joint 
are  the  trochlea  and  the  trochlear  notch.  The  trochlea 
(Figure  18a)  is  a  grooved  surface  much  like  the  circumference 
of  a  pully,  which  covers  the  anterior,  inferior  and  posterior 
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Table  1 


DEGREE  OF  POLYNOMIAL  (n) ,  STANDARD  DEVIATION  (o) 
AND  PERCENTAGE  ACCURACY  (%R)  OF  THE  FIT  FOR 
THE  ARTICULATING  SURFACES  OF  ELBOW,  HIP, 
KNEE  AND  ANKLE  JOINTS. 


JOINT 

ARTICULATING 

SURFACE 

n 

a  (cm) 

°sR 

ELBOW 

Trochlea 

4 

0.09 

94.4 

Trochlear  notch 

4 

0.07 

94 . 5 

HIP 

Head  of  femur 

4 

0.08 

98.2 

Acetabulum 

4 

0.15 

99.1 

KNEE 

Tibia  lateral 

4 

0.03 

95.8 

Tibia  medial 

4 

0.05 

94.6 

Femur  lateral 

4 

0.04 

99.3 

Femur  medial 

4 

0.07 

98.7 

ANKLE 

Talus  (trochlear) 

4 

0.03 

98.5 

Medial  malleolus 

4 

0.04 

95.7 

surfaces  of  the  condyle  of  the  humerus.  It  is  separated  from 
the  capitulum  on  its  lateral  side  by  a  faint  groove,  but  its 
medial  margin  is  salient  and  projects  downward  beyond  the 
rest  of  the  bone.  The  trochlea  articulates  with  the  trochlear 
notch  of  the  ulna. 

The  trochlear  notch  (Figure  18b)  is  formed  by  the  anterior 
surface  of  the  olecranon  and  the  superior  surface  of  the 
coronoid  process.  The  base  is  constricted  at  the  junction 
between  these  two  areas  and  they  may  be  separated  completely 
by  a  narrow,  roughened  strip.  A  smooth  ridge  which  corresponds 
to  the  groove  of  the  trochlea,  divides  the  notch  into  a  lar.ger, 


medial  portion  and  a  smaller,  lateral  part.  The  medial  part 
conforms  to  the  large  flange  of  the  trochlea  of  the  humerus. 

Coordinate  System 

The  origin  of  the  (x,y,z)  coordinate  system  is  placed  at 
the  approximate  geometric  center  of  the  humerus,  with  the 
x-axis  directed  along  the  posterior-anterior  direction  and 
the  y-axis  Coinciding  with  the  humerus  longitudinal  axis 
(Figure  18a).  The  origin  of  the  coordinate  system  (x',y',z') 
coincides  approximately  with  the  geometric  center  of  trochlear 
notch,  with  the  x'-axis  directed  along  the  anterior-posterior 
direction  and  the  y'-axis  being  directed  along  the  longitudi¬ 
nal  axis  of  the  ulna  (Figure  18b)  .  The  locations  of  the 
origins  of  these  coordinate  systems  are  15.5  cm  and  -0.8  cm 
from  the  intersection  points  of  the  y  and  x'  axes  with  the 
articular  surfaces,  respectively. 

The  articulating  surfaces  are  digitized  in  their 
respective  coordinate  systems  following  the  technique  described 

above.  The  coefficients  a  of  equation  (1)  are  summarized  in 

pq 

Table  2.  Note  that  to  avoid  multi-value  function  difficulties, 
the  surface  equation  for  the  trochlear  notch  is  presented  as: 

x 1  =  f (y ' , z ' )  -  l  a  y,pz'q  (2) 

p=0  q=0  Pq 

Proper  consideration  of  this  variable  change  must  be  made  in 
future  analysis  and  modeling  of  the  elbow  joint. 

ARTICULATING  SURFACES  OF  THE  HIP  JOINT 
Anatomical  Description 

The  two  articulating  surfaces  of  the  hip  joint  are  the 
head  of  the  femur  and  the  acetabulum  of  the  pelvis.  The  head 
of  the  femur  is  rather  more  than  half  of  a  sphere  (Figure  19a). 
It  is  directed  upward,  medially  and  slightly  forward  to 
articulate  with  the  acetabulum.  Its  surface  is  smooth  with  a 
small  fovea  or  roughened  pit  slightly  below  and  behind  its 
center. 
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Table  2 


COEFFICIENTS  OF  THE  EQUATIONS  OF  THE 
ELBOW  JOINT  ARTICULATING  SURFACES. 


COEFFICIENTS 

apq 

(cm) 

ARTICULATING  SURFACE 

Trochlea 

Trochlear  Notch 

a00 

15.5619 

-0.8296 

a01* 

0.5426 

0.1613 

a02 

0.9899 

0.0410 

a03 

-0.3204 

-0.2188 

a04 

-0.3807 

0.1936 

aio* 

0.5668 

0.4230 

all 

-0.3387 

-0.2996 

a12 

0.1875 

-0.1772 

aI3 

0.2360 

0.3742 

a20 

-0.4083 

-0.5579 

a21 

-0.0320 

-0.0439 

a22 

-0.2297 

-0.7558 

a30 

-0.0243 

-0.2487 

a31 

0.1307 

-0.0776 

a40 

-0.0891 

0.1855 

*unitless 
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Figure  19.  Anterior  and  front  aspects  of  the  femur 
and  pelvis. 


-43- 


The  acetabulum  is  an  approximately  hemispherical  cavity 
on  the  lateral  asoect  of  the  inominate  bone,  about  its  center, 
and  is  directed  laterally,  downward  and  forward.  It  is 
surrounded  by  an  irregular  projecting  margin  which  is  defined 
inferiorly;  this  gap  is  the  acetabular  notch.  The  floor  of 
the  cavity  is  roughened  and  non-articular .  The  sides  of  the 
cup  present  an  articular  lunate  surface  which  is  widest 
superiorly.  In  this  situation,  the  weight  of  the  trunk  is 
transmitted  to  the  femur  in  the  erect  attitude.  This  horse¬ 
shoe  shaped  strip  is  covered  with  articular  cartilage  (Gray 
[1973])  and  provides  the  surface  on  which  the  head  of  the 
femur  rides  within  the  hip  joint. 

Coordinate  System 

The  origin  of  the  (x,y,z)  coordinate  system  is  placed  at 
the  geometric  center  of  the  acetabulum  floor  (acetabular 
fossa) ,  with  the  x-axis  directed  along  the  posterior-anterior 
direction  and  the  y-axis  directed  along  the  superior- inferior 
direction  (Figure  19b).  The  origin  of  the  (x’,y',z') 
coordinate  system  coincides  with  the  approximate  center  of 
mass  of  the  femur,  with  the  x'-axis  directed  along  the 
posterior-anterior  direction  and  the  y'-axis  being  directed 
along  the  longitudinal  axis  of  the  femur  (Figure  19a) .  The 
location  of  the  origins  of  these  coordinate  systems,  as  shown 
in  Figure  19,  are  -1.5  cm  and  19.5  cm  from  the  intersection 
points  of  the  y  and  y'  axes  with  the  articular  surfaces, 
respectively. 

The  articulating  surfaces  of  the  hip  joint  are  digitized 
in  their  respective  coordinate  systems  as  described  previously 
and  the  results  for  the  coefficients  ap^  of  equation  (1)  are 
summarized  in  Table  3. 


ARTICULATING  SURFACES  OF  THE  KNEE  JOINT 
Anatomical  Description 

The  tibial  and  femoral  condyles  are  the  major  articulating 
surfaces  of  the  knee  joint  (Figure  20).  The  upper  end  of  the 


Table  3 

COEFFICIENTS  OF  THE  EQUATIONS  OF  THE 
HIP  JOINT  ARTICULATING  SURFACES. 


COEFFICIENTS 

a 

ARTICULATING 

SURFACE 

Pm 

(cm) 

Head  of  Femur 

Acetabulum 

a00 

-19.5628 

-1.5131 

a01* 

-3.9686 

-1.2324 

a02 

1.7498 

-0.0416 

a03 

-0.3864 

0.5052 

a04 

0.0352 

0.1025 

a10* 

0.2045 

4.5526 

all 

0.2260 

10.5723 

a12 

-0.1281 

4.9788 

a13 

0.0174 

0.6813 

a20 

0.9641 

-0.0956 

a21 

-0.4715 

-0.4650 

a2  2 

0.0800 

0.0034 

a30 

-0.0382 

1.5829 

a31 

0.0146 

0.5512 

a40 

0.0235 

-0.1303 

Figure  20.  Posterior  aspect  of  the  femur  and  anterior 
view  of  the  tibia. 
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tibia  is  expanded,  especially  in  its  transverse  axis,  providing 
an  adequate  bearing  surface  for  the  body  weight  transmitted 
through  the  lower  end  of  the  femur.  It  comprises  of  two  prom¬ 
inent  masses,  the  medial  and  lateral  condyles,  and  a  smaller 
projection,  the  tuberosity  of  the  tibia.  The  condyles  project 
backwards  a  little,  so  as  to  overhang  the  upper  part  of  the 
posterior  surface  of  the  shaft.  Superiorly,  each  is  covered 
with  an  articular  surface,  the  two  being  separated  by  an 
irregularly  roughened  intercondylar  area.  They  form  visible 
and  palpable  landmarks  at  the  side  of  the  ligamentum  patella, 
the  lateral  condyle  being  the  more  prominent. 

The  medial  condyle  (Figure  20b)  is  the  larger  but  does 
not  overhang  so  much  as  the  lateral  condyle.  Its  upper 
articular  surface,  oval  in  outline,  is  concave  in  all 
diameters,  and  its  lateral  border  projects  upwards,  deepening 
the  concavity  and  covering  an  elevation,  the  medial  intercon¬ 
dylar  tubercle.  The  posterior  surface  of  the  condyle  is 
marked,  immediately  below  the  articular  margin,  by  a 
horizontal,  roughened  groove.  Its  medial  and  anterior 
surfaces  form  a  rough  strip,  separated  from  the  medial  sur¬ 
face  of  the  shaft  by  an  inconspicuous  ridge. 

The  lateral  condyle  (Figure  20b)  overhangs  the  shaft, 
especially  at  its  posterolateral  part,  which  bears  on  its 
inferior  surface  a  small  circular  facet  for  articulation  with 
the  upper  end  of  the  fibula.  The  upper  surface  is  covered 
with  an  articular  surface  for  the  lateral  condyle  of  the 
femur.  Nearly  circular  in  outline,  it  is  slightly  hollowed 
in  its  central  part,  and  its  medial  border  extends  upwards  to 
cover  an  elevation,  termed  the  lateral  intercondylar  tubercle. 
The  posterior,  lateral  and  anterior  surfaces  of  the  condyle 
are  rough. 

The  lower  end  of  the  femur  is  widely  expanded  and  thu? 
provides  a  good  bearing  surface  for  the  transmission  of  the 
weight  of  the  body  to  the  top  of  the  tibia.  It  consists  of 


two  prominent  masses  of  bone,  the  condyles  (Figure  20a), 
which  are  partially  covered  by  a  large  articular  surface. 
Anteriorly,  the  two  condyles  are  united  and  are  continuous 
with  the  front  of  the  shaft;  posteriorly,  they  are  separated 
by  a  deep  gap,  the  intercondylar  fossa  (intercondylar  notch), 
and  they  project  backwards  considerably  beyond  the  plane  of 
the  popliteal  surface. 

The  lateral  condyle  is  flattened  on  its  lateral  surface 
and  is  not  so  prominent  as  the  medial  condyle,  but  it  is 
stouter  and  stronger,  for  it  is  placed  more  directly  in  line 
with  the  shaft  and  probably  takes  a  greater  share  in  the 
transmission  of  the  weight  to  the  tibia.  The  most  prominent 
point  on  its  lateral  aspect  is  termed  the  lateral  epicondyle 
(Figure  20a)  ,  and  the  whole  of  this  surface  can  be  felt 
through  the  skin. 

The  medial  condyle  possesses  a  bulging,  convex  medial 
aspect,  which  can  be  palpated  without  difficulty.  Its  upper¬ 
most  part  is  marked  by  a  small  projection,  termed  adductor 
tubercle  because  it  gives  insertion  to  the  tendon  of  the 
adductor  magnus .  The  most  prominent  point  on  the  medial  sur¬ 
face  of  the  condyle  is  below  and  a  little  in  front  of  the 
adductor  tubercle  and  is  termed  the  medial  epicondyle  (Figure 
20a) .  The  lateral  surface  of  the  condyle  is  the  roughened 
medial  wall  of  the  intercondylar  fossa. 

Coordinate  Systems 

The  origin  of  the  (x,y,z)  coordinate  system  is  placed  at 
the  approximate  geometric  center  of  the  femur  with  the  x-axis 
directed  along  the  posterior-anterior  direction  and  the  y-axis 
coinciding  with  the  femural  longitudinal  axis  (Figure  20a). 

The  origin  of  the  coordinate  system  (x',y',z')  coincides  with 
the  approximate  center  of  mass  of  the  tibia,  with  the  x'-axis 
directed  along  the  anterior-posterior  direction  and  the 
y'-axis  being  directed  along  the  longitudinal  axis  of  the 
tibia  (Figure  20b) .  The  location  of  the  origins  of  these 
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coordinate  systems,  shown  in  Figure  20,  are  19.5  and  -20.6  cm 
from  the  intersection  points  of  y  and  y*  axes  with  the  lateral 
articular  surfaces. 

The  articulating  surfaces  of  the  knee  joint  are  digitized 
in  their  respective  coordinate  systems  according  to  the 
technique  described  previously  and  the  coefficients,  a  ,  of 

r4 

equation  (1)  for  the  tibial  and  femoral  articulating  surfaces 
of  the  knee  joint  are  summarized  in  Tables  4  and  5, 
respectively. 


ARTICULATING  SURFACES  OF  THE  ANKLE  JOINT 
Anatomical  Description 

The  major  articulating  surfaces  of  the  ankle  joint  are 
the  trochlear  surface  of  the  talus  and  the  medial  malleolus 
(Figure  21) . 

The  body  of  the  talus  is  cuboidal  in  shape.  Its  dorsal 
surface  is  covered  by  the  trochlear  articular  surface,  which 
articulates  with  the  lower  end  of  the  tibia  at  the  ankle 
joint.  It  is  convex  from  back  to  front  and  gently  concave 
from  side  to  side,  and  it  is  widest  anteriorly  (Figure  21a). 
The  medial  surface  of  the  talus  is  covered  in  its  upper  part 
by  a  comma-shaped  articular  facet  which  is  deeper  in  front 
than  behind  and  articulates  with  the  medial  malleolus. 

The  medial  malleolus  is  a  short  but  stout  process.  Its 
lateral  surface  is  smooth  and  occupied  by  a  comma-shaped 
articular  facet,  which  articulates  with  the  medial  side  of 
the  talus  (Figure  21b).  Its  anterior  surface  is  rough,  and 
its  posterior  surface  bears  the  lower  end  of  the  groove  that 
marks  the  posterior  surface  of  the  lower  end  of  the  bone. 

The  lower  border  of  the  malleolus  is  pointed  anteriorly  and 
depressed  posteriorly. 

Coordinate  Systems 

The  origin  of  (x,y,z)  coordinate  system  is  placed  at  the 
approximate  geometric  center  of  the  tibia,  with  the  x-axis 


Table  4 


COEFFICIENTS  OF  THE  EQUATIONS  OF  THE  TIBIAL 
ARTICULATING  SURFACES  OF  THE  KNEE  JOINT. 


COEFFICIENTS 

a 

pq 

(cm) 

ARTICULATING  SURFACE 

Tibia  Lateral 

Tibia  Medial 

a00 

-20.6964 

-18.7135 

a0l* 

-1.6128 

-1.6587 

a02 

-0.3786 

1.9298 

a03 

0.0215 

-0.7677 

a04 

0.0117 

0.1001 

‘  a10* 

0.7539 

-0.7837 

S11 

1.0017 

0.9363 

a12 

0.3751 

-0.4480 

a13 

0.0421 

0.0677 

a20 

0.2212 

-0.3921 

a21 

0.1889 

0.1772 

a2  2 

0.0335 

-0.0123 

a30 

0.0885 

-0.2331 

a31 

0.0367 

0.0670 

a40 

0.0373 

-0.0369 

*unitless 
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Table  5 


COEFFICIENTS  OF  THE  EQUATIONS  OF  THE  FEMORAL 
ARTICULATING  SURFACES  OF  THE  KNEE  JOINT. 


COEFFICIENTS 

apq 

(cm) 


40 


*unitless 


ARTICULATING  SURFACE 

Femur  Lateral 

Femur  Medial 

19.5520 

18.5270 

-1.9895 

8.9256 

0.6298 

-7.2093 

0.8031 

2.4873 

-0.0055 

-0.3193 

-1.4614 

0.1211 

-2.5323 

-0.5608 

-1.3528 

0.1310 

0.2354 

-0.0500 

-0.1912 

-0.5928 

0.0620 

-0.0411 

-0.4573 

-0.0557 

-0.1305 

-0.2447 

0 . 2058 

-0.0332 

-0.0544 

-0.0376 
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Figure  21.  Medial  and  posterior  aspects  of  the  talus 
and  malleolus. 
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directed  along  the  posterior-anterior  direction  and  the  y-axis 
coinciding  with  the  longitudinal  axis  of  tibia  (Figure  21b) . 
The  origin  of  the  coordinate  system  (x',y',z')  coincides  with 
the  approximate  center  of  mass  of  the  talus,  with  the  x'-axis 
directed  along  the  posterior-anterior  direction  and  the  y ' - 
axis  being  directed  along  the  inferior- superior  direction  of 
the  talus  (Figure  21a) .  The  location  of  the  origins  of  these 
coordinate  systems,  shown  in  Figure  21,  are  20.2  cm  and  1.3  cm 
from  the  intersection  points  of  the  y  and  y*  axes  with  the 
articular  surfaces,  respectively. 

The  articulating  surfaces  of  the  ankle  joint  are 
digitized  in  their  respective  coordinate  systems  according  to 
the  technique  described  previously  and  the  coefficients,  a  , 

r4 

of  equation  (1)  are  summarized  in  Table  6. 

MATHEMATICAL  REPRESENTATION  OF  LIGAMENTS 

Structural  integrity  of  the  articulating  joints  is 
maintained  by  capsular  ligaments  and  both  extra-  and  intra- 
articular  ligaments.  Capsular  ligaments  are  formed  by  thick¬ 
ening  of  the  capsule  walls  where  functional  demands  are 
greatest.  As  the  names  imply,  extra-  and  intra-articular 
ligaments  at  the  joints  reside  external  to  and  internal  to  the 
joint  capsule,  respectively.  Extra-articular  ligaments  have 
several  shapes,  e.g.  cord-like  or  flat  depending  on  their 
locations  and  functions.  These  types  of  ligaments  appear 
abundantly  at  the  articulating  joints.  However,  only  the 
shoulder,  hip  and  knee  joints  contain  intra-articular 
ligaments.  For  example,  the  cruciate  ligaments  at  the  knee 
joint  are  probably  the  most  well  known  intra-articular 
ligaments.  Further  information  about  the  structure  and 
mechanics  of  the  joint  can  be  found  in  Barnett,  Davies,  and 
Mac  Conaill  [1949]. 

v 

Ligament  configurations  are  largely  dependent  on 
arrangement  of  the  joint  and  its  articulation,  the  direction 
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COEFFICIENTS  OF  THE  EQUATIONS  OF  THE  ARTICULATING 
SURFACES  OF  THE  ANKLE  JOINT. 


COEFFICIENTS 

apq 

(cm) 

ARTICULATING  SURFACE 

Talus  (Trochlear) 

Medial  Malleolus 

P 

o 

o 

-1.4684 

20.2305 

a01* 

0.3722 

-0.0445 

a02 

-0.1347 

-0.1872 

a03 

-0.0227 

0.0010 

a04 

0.0549 

0.0601 

aio* 

0.1428 

-0.0689 

all 

0.0018 

0.1748 

a12 

-0.1114 

-0.0195 

a13 

-0.0454 

-0.0725 

a20 

0.3446 

0.1118 

a21 

-0.0123 

0.0319 

a22 

-0.0115 

0.0757 

a30 

0.0516 

0.0369 

a31 

0.0154 

-0.0232 

a40 

-0.0072 

0.0316 

*unitless 


of  tendons,  and  the  location  of  the  organs.  Thus,  in 
different  joints,  according  to  their  structure,  ligaments  may 
be  more  or  less  lax,  or  tight,  or  capable  of  more  or  less 
movement.  Because  of  the  difference  in  the  direction  of  the 
tendons  or  the  placement  of  the  membranes  and  their  site  and 
magnitude,  ligaments  vary  in  their  configurations,  circumvolu¬ 
tion  and  extension.  On  account  of  all  these  factors,  some 
ligaments  are  twisted  like  cords,  some  are  united  in  fibrous 
bonds  and  others  flattened  into  membranous  forms.  This 
variation  in  form  is  indicated  in  their  corresponding  names 
according  to  size:  major,  minor,  maximus;  according  to 
external  form:  large,  thick,  thin;  according  to  configuration: 
long,  wide,  round,  triangular,  quadrate,  circular;  according 
to  their  positions:  straight,  transverse,  oblique,  horizontal, 
perpendicular,  superficial,  sublime,  deep,  lateral,  right, 
left,  anterior,  posterior,  superior  and  inferior;  according 
to  their  insertion:  interclavicular ,  brachioradial ,  etc. 

(Gray  [1973]). 

In  this  section,  general  aspects  and  material 
characteristics  of  soft  tissues,  in  general,  and  ligaments  in 
particular,  are  studied.  Ligaments  are  modeled  as  non-linear 
elastic  springs  and  their  constitutive  equation  and  stiffness 
values  are  presented.  Attachment  sites  of  the  elbow,  hip, 
knee  and  ankle  joint  ligaments  are  determined  and  summarized. 

GENERAL  CHARACTERISTICS 

Ligaments,  capsule  and  other  connective  tissues  such  as 
tendons,  skin  and  blood  vessels,  consist  mainly  of  collagen 
and  elastin  fibers  embedded  in  a  mucopolysaccharide  intercel¬ 
lular  ground  substance  (Wismans  [1980]).  Geometrical  arrange¬ 
ments  and  the  relative  amount  of  each  of  the  components  of 
the  fibers  vary  from  tissue  to  tissue.  Usually,  in  ligaments 
and  tendons,  collagen  fibers  are  oriented  in  the  direction  of 
the  transmitted  force  while  the  elastin  fibers  form  a  dis¬ 
ordered  network.  Crisp  [1972]  reported  that  under  no  external 
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loading  of  the  tissue,  the  collagen  bundles  are  coiled.  Due 
to  the  difficulties  involved  in  the  dissection  of  the  compo¬ 
nents  of  biological  structures,  the  mechanical  properties  of 
the  collagen  and  elastin  fibers  are  usually  determined  for  a 
tissue  in  which  one  of  the  components  is  predominent. 

Material  properties  of  collagen  are  often  determined  by  the 
use  of  tendons.  Elliot  [1965]  reported  that  in  tendons  of 
the  human  body,  about  75%  of  the  dry  weight  is  collagen  and 
just  21  is  elastin. 

Physical  behavior  and  the  stress  in  a  biological  tissue 
not  only  depends  on  strain  but  also  on  strain  history.  This 
behavior  is  evident  in  phenomena  such  as  stress  relaxation, 
creep,  hysteresis  and  dependence  of  the  elastic  moduli  on  the 
strain  rate  and  temperature.  A  number  of  mathematical 
descriptions  characterizing  this  behavior  are  proposed  in  the 
literature.  ese  descriptions  may  be  divided  into  two  groups. 

In  the  first  group,  the  measured  microscopic  response  of 
the  tissue  is  characterized  by  a  continuous,  nonlinear 
equation,  as  in  the  quasi-linear  visco-elastic  law  of  Fung 
[1972].  The  most  noticeable  feature  of  the  mechanical 
behavior  of  biological  tissues  is  that  measurable  stresses 
develop  only  after  the  specimen  has  been  stretched  consider¬ 
ably  from  its  original  or  relaxed  length.  In  such  an 
extension,  the  stress-strain  law  becomes  highly  nonlinear  and 
classical  theory  of  elasticity,  which  is  restricted  to  linear 
stress-strain  relations  and  small  strains,  is  not  applicable. 
The  nonlinear  stress  -  strain  relations  developed  in  the 
analysis  of  finite  homogeneous  deformations  of  elastomers 
have  been  studied  by  Mooney  [1940],  Rivlin  [1948]  and  on  the 
basis  of  strain  energy  function  by  Green  and  Adkins  [I960]. 

Fung  [1967]  has  shown  that  the  elastic  properties  of  mesentery 
are  completely  different  from  those  of  vulcanized  rubber.  He 
has  concluded  that  the  stress  -  strain  relation  in  the  one¬ 
dimensional  case  should  be  exponential  in  the  stretch  mode. 


A  generalization  of  Fung's  result  to  three-dimensional 
problems  was  given  by  Gou  [1970],  who  introduced  a  strain 
energy  density  which  is  an  exponential  function  of  the  strain 
invariant.  Various  potential  functions  describing  the  large 
deformation  of  biological  tissue  and  its  response  to  various 
forces  are  presented  in  the  literature.  These  include  the 
works  by  Blatz,  Chu  and  Wayland  [1969],  Lee,  Frasher  and  Fung 
[1967],  Hildebrandt,  Fukaya  and  Martin  [1969],  Veronda  and 
Westman  [1970],  Simon,  et  al.,  [1970],  and  Demiray  [1972]. 

In  general,  however,  these  functions  have  been  developed  for 
a  specific  loading  pattern  and  assumed  material  isotropy.  A 
more  general  function  permitting  the  study  of  a  number  of 
types  of  loads  and  interaction  between  combined  loads  is  given 
by  Snyder  [1972] . 

In  this  first  group,  there  are  also  represented  models 
in  which  the  tissue  response  is  described  by  a  mechanical 
analogy,  consisting  of  a  number  of  spring,  dashpot  and  dry 
frictional  elements.  For  example,  nonlinear  visco-elastic 
behavior  of  collagenous  tissue  has  been  simulated  by  Frisen, 
et  al.,  [1969],  using  a  Kelvin  model  and  a  number  of  nonlinear 
springs . 

In  the  second  group  the  mathematical  description  is  based 
on  an  idealization  of  the  microstructures  and  on  the  mechani¬ 
cal  properties  of  the  constituent  materials.  Based  on  the 
knowledge  that  the  connective  tissues  are  composed  of  fibrous 
and  amorphous  materials  they  may  be  treated  as  fiber-reinforced 
materials.  The  initial  part  of  the  loading  phase  is  a  geomet¬ 
rical  rearrangement  of  the  microstructural  network  due  to 
uncoiling  of  the  coiled  fibers  of  collagen.  Linear  constitu¬ 
tive  equations  are  assumed  for  the  second  part  of  the  loading 
in  which  the  collagen  fibers  are  elongated.  The  models  by 
Comninou  and  Yannas  [1976]  ,  and  Drouin  [1980]  on  the  one¬ 
dimensional  stress  field  studies  of  elastic  behavior  of 


-57- 


collagenous  tissues  and  two  material  composite  prosthesis, 
respectively,  fall  in  this  group. 

In  this  study  simulation  of  joint  ligaments  is 
accomplished  by  a  mathematical  description  based  on  an 
approximation  of  experimental  data  reported  in  the  literature. 
Ligaments  will  be  represented  by  nonlinear  elastic  springs  as 
the  experiments  were  limited  to  a  one-dimensional  stress 
field  and  information  is  available  only  on  the  elastic 
behavior.  Due  to  the  microstructure  and  the  orientation  of 
the  collagen  fibers,  representation  of  ligaments  by  spiings 
would  seem  to  be  a  realistic  approach.  The  general  form  of 
the  constitutive  equation  and  the  values  of  the  parameters 
characterizing  the  constitutive  equations  of  the  different 
springs  are  discussed  below. 

CONSTITUTIVE  EQUATION 

A  constitutive  equation  representing  ligamentous  behavior 
is  based  on  available  data  in  the  literature.  Data  concerned 
with  the  knee  joint  ligaments  is  considered  since  the  general 
theoretical  analysis  developed  later  in  this  report  will  be 
applied  to  the  knee  joint. 

Brantigan  and  Voshell  [1941]  presented  a  review  of  the 
conflicting  theories  on  the  function  of  knee  ligaments  prior 
to  1940  and  reported  the  results  of  study  on  approximately  100 
knees.  Since  that  time  many  investigators  have  discussed  the 
function  of  various  ligamentous  structures  (Hallen  and  Lindhal 
[1965]  and  [1966];  Hughston  and  Eilers  [1973];  Kennedy  and 
Grainger  [1967];  Kennedy  and  Fowler  [1971];  Kennedy,  Weinberg 
and  Wilson  [1974];  Robinson  and  Romero  [1968];  Slocum  and 
Larson  [1968];  Warren,  Marshall  and  Girgis  [1974];  Girgis,  et 
al.,  [1975];  Trent  and  Walker,  [1975];  Piziali,  Rastegar  and 
Nagel  [1977];  Piziali,  et  al.,  [1980];  and  Seering,  et  al., 
[1980]),  and  the  length  of  primary  structures  as  a  function 
of  knee  flexion  (Edwards,  Lafferty  and  Lange  [1969];  Wang, 
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Walker  and  Wolf  [1973];  and  Crowninshield,  Pope  and  Johnson 
[1976].  In  addition,  pure  torsional  rotation  with  and  without 
contact  pressure  has  been  measured  by  Wang  and  Walker  [1974]. 

Force -deflexion  characteristics  of  several  types  of 
ligaments  under  uniaxial  tension  has  been  reported  by  Trent, 
Walker  and  Wolf  [1976],  Kennedy,  et  al.,  [1976],  Noyes  and 
Grood  [1976]  and  Dorlot,  et  al.,  [1980].  Usually  four 
characteristic  regions  of  the  force-deflexion  curve  of  collag¬ 
enous  tissue  are  identified.  The  location  of  these  regions 
(Figure  22)  can  be  explained  in  terms  of  microarchitecture 
(Crisp  [1972];  Wismans  [1980]).  The  stiffness  of  the  tissue, 
defined  by  the  slope  of  the  load-deflexion  curve,  is  rather 
slight  in  region  1.  This  initial  region  is  considered  mainly 
to  correspond  to  the  geometrical  rearrangement  of  the  micro- 
structural  network  (uncoiling  of  the  coiled  collagen  fibers). 
Therefore,  the  stiffness  is  determined  mainly  by  the  stiffness 
of  the  elastin  network.  The  stiffness  increases  in  region  2 
as  some  of  the  fibers  become  aligned.  All  collagenous  fibers 
are  assumed  to  be  fully  uncoiled  at  the  end  of  this  region. 

As  the  force  on  the  tissue  steadily  increases,  the  collagen 
fibers  themselves  elongate.  In  region  3,  the  stiffness  is 
reported  to  correspond  mainly  to  the  stiffness  of  the  collagen 
fibers  and  is  found  to  be  almost  constant.  Finally  in  region 
4,  disruption  of  some  collagen  fibers  is  observed,  followed 
by  a  complete  failure  of  the  tissue  itself. 

The  force-elongation  curve  of  Figure  22  can  be 
represented  by  a  quadratic  equation.  Haut  and  Little  [1972] 
carried  out  a  number  of  tension  tests  on  rat-tail  collagen 
bundles  and  reported  that  at  low  strains,  the  elastic  behavior 
of  this  tissue  could  be  described  by  a  quadratic  stress-strain 
equation.  Crowninshield,  Pope  and  Johnson  [1976]  tested  human 
medial  collateral  ligaments  and  indicated  that  a  quadratic 
stress-strain  function  is  a  good  approximation  for  the  elastic 


-59- 


STRAIN!,  * 
LENGTH. I 


disruption  of 
coiiogen  fibers 


Figure  22.  Force-elongation  curve  for  collagenous  tissues. 

behavior  of  such  tissue.  In  this  study  the  following  force- 
elongation  relationship  is  assumed  for  each  ligament: 


fj|  ‘  Vlj  -  V 


for  Lj  > 


in  which  kj  is  the  spring  constant,  Lj  and  l .  are,  respectively 
the  current  and  initial  lengths  of  the  ligament,  j.  It  is 
assumed  that  the  ligaments  cannot  carry  any  compressive  force; 
accordingly: 


fjl  -  0 


for  Lj  <  tj 


The  direction  of  the  force,  Fj  ,  exerted  by  a  spring  on  the 
articulating  body  segment  coincides  with  the  direction  of  the 
line  segment  through  the  origin  and  insertion  points  of  that 
spring.  The  length,  l . ,  of  spring  j  is  equal  to  the  distance 
between  its  insertion  and  origin  points. 


LIGAMENT  STIFFNESS 

Several  studies  concerned  with  the  tension  tests  of  human 
knee  joint  ligaments  have  been  reported  in  the  literature. 
Among  these  studies,  the  work  of  Noyes  and  Grood  [1976], 

Trent,  Walker  and  Wolf  [1976],  and  Kennedy,  et  al.,  [1976] 
will  be  considered  in  this  section. 

Twenty-one  anterior  cruciate  ligaments  have  been  tested 
by  Noyes  and  Grood  [1976].  The  intact  knee  joints  were  kept 
in  a  frozen  condition  for  about  four  weeks.  They  were  thawed 
at  room  temperature  a  day  before  testing.  The  specimens  were 
then  mounted  in  the  testing  machine  at  45°  of  flexion  and,  at 
an  elongation  rate  of  25  mm/sec.,  they  were  tested  to  failure. 
The  stiffness,  k,  the  strain  e £,  and  the  corresponding  force, 
F£,  were  determined.  The  mean  values  of  these  quantities  are 
summarized  in  Table  7.  Large  deviations  between  the  specimens 
are  evident  and  the  younger  ligaments  are  much  stiffer  than 
older  ones. 

The  cruciates  and  the  collateral  ligaments  of  six  fresh 
specimens  together  with  a  piece  of  their  bony  attachments  were 
tested  by  Trent,  Walker  and  Wolf  [1976].  The  specimens  were 
loaded  to  failure  in  a  normal  saline  solution  at  37°C  and  at 
an  elongation  rate  of  0.8  mm/sec.  Table  7  shows  the  mean- 
values  for  the  stiffness,  k.  The  values  for  the  anterior 
cruciate  is  close  to  the  mean  values  reported  by  Noyes  and 
Grood  [1976]  for  the  older  specimens. 

Kennedy,  et  al.,  [1976]  tested  twenty  anterior,  posterior 
and  collateral  ligaments  at  two  different  strain  rates  of  2 
and  8  mm/second.  The  specimens  were  tested  some  fourteen 
hours  after  death  and  were  kept  in  isotonic  saline.  The 
ligaments  were  tested  without  their  bony  attachments.  Each 
ligament  was  fixed  in  an  upper  clamp  initially  and  allowed  to 
hang  freely,  seeking  its  own  orientation.  A  lower  clamp  was 
then  applied.  Special  care  was  made  to  avoid  twisting  the 
ligaments.  The  mean  values  of  the  strain,  e£,  and 
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Table 


corresponding  force,  F£,  are  given  m  Table  7.  No  stiffness 
data  or  age  of  specimens  were  reported  in  this  study.  The 
variation  in  the  results  reported  by  Kennedy,  et  al.,  [1976] 
and  those  of  Noyes  and  Grood  [1976]  may  be  explained  partly 
by  possible  slipping  of  the  specimens  between  the  clamps  and 
differences  in  elongation  rates. 

In  this  study,  the  values  of  the  stiffness,  kj ,  for 
different  springs  is  based  on  the  mean  data  reported  by  Trent, 
Walker  and  Wolf  [1976].  For  the  stiffness  of  the  medial 
collateral  ligament  a  correction  reported  by  Wisman  [1980]  is 
used  in  order  to  account  for  the  oblique  and  the  deep  medial 
collateral  ligament,  since  these  parts  were  not  included  in 
the  specimens  tested  by  Trent,  Walker  and  Wolf  [1976].  This 
correction  is  based  on  data  presented  by  crowninshield ,  Pope 
and  Johnson  [1976].  Stiffness  values  used  in  this  study  for 
the  lateral  collateral  (LC)  ,  the  medial  collateral  (MC) ,  the 
anterior  cruciate  (AC),  and  the  posterior  cruciate  (PC) 
ligaments  are  summarized  in  Table  8. 


Table  8 

LIGAMENT  STIFFNESS  FOR  THE  LATERAL 
COLLATERAL  (LC)  ,  MEDIAL  COLLATERAL 
(MC),  ANTERIOR  CRUCIATE  (AC)  AND 
POSTERIOR  CRUCIATE  (PC)  LIGAMENTS. 


Ligament 
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INSERTIONS  AND  ORIGINS 

Accurate  determination  of  attachment  sites  or  the  so- 
called  insertions  and  origins,  of  ligaments  are  of  considerable 
importance  in  the  mathematical  modeling  of  the  articulating 
joint.  Very  little  quantitative  information  is  available  on 
this  subject.  Attachment  sites  and  the  lengths  of  the  knee 
joint  ligaments  have  been  measured  in  vivo  and  in  vitro  and 
reported  by  Crowninshie Id ,  Pope,  and  Johnson  [1976]. 

In  this  study,  by  close  physical  and  anatomical  study  of 
the  elbow,  hip,  knee,  and  ankle  joints  (Grant  [1962];  Gray 
[1973]),  the  attachment  sites  of  their  respective  ligaments 
have  been  measured  and  summarized  in  Table  9.  For  each  joint, 
insertion  points  are  measured  with  respect  to  its  (x',y',z') 
coordinate  system  and  origins  with  respect  to  its  (x,y,z) 
coordinate  system.  Note  that  attachment  sites  of  only  those 
ligaments  which  significantly  contribute  to  the  integrity  and 
functional  aspects  of  the  joints  have  been  measured,  and  due 
to  the  expanded  shape  of  some  of  the  ligaments,  slight  varia¬ 
tion  may  be  inherent  in  these  measurements. 

TWO-DIMENSIONAL  DYNAMIC  FORMULATION 
OF  A  TWO -BODY  SEGMENTED  JOINT 

The  mathematical  descriptions  of  the  articular  surfaces, 
ligaments  and  capsule  thus  far  presented  can  now  be  integrated 
into  a  two-dimensional  mathematical  formulation  of  a  general 
two-body,  segmented  articulating  joint. 

For  the  purpose  of  studying  the  joint  motion,  one  segment 
is  assumed  to  be  fixed  while  the  other  segment  is  executing  a 
relative  motion.  The  coordinate  systems  (x,y)  and  (x',y')  are 
attached  to  the  fixed  and  the  moving  body  segments, 
respectively,  and  their  relative  position  and  angular  orienta¬ 
tion  will  be  discussed.  Next  the  contact  conditions  between 
the  two  articulating  surfaces  are  presented,  and  the  descrip¬ 
tion  of  ligaments,  contact  and  applied  external  forces  and 
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INSERTIONS  AND  ORIGINS  FOR  THE  LIGAMENTS  OF  THE 
ELBOW,  HIP,  KNEE  AND  ANKLE  JOINTS 
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moments,  along  with  the  governing  dynamic  equations  of  motion 
are  presented.  A  numerical  procedure  employing  the  Newmark 
method  of  differential  approxime: tion  and  Newton-Raphson  itera 
tion  process  is  suggested  for  the  solution  of  these  coupled, 
nonlinear  algebraic  and  differential  equations. 

CHARACTERIZATION  OF  THE  RELATIVE  POSITIONS 

A  joint  connects  two  segments  of  a  body  which  are 
designated  as  segments  1  and  2  and,  for  illustrative  purposes 
will  be  represented  as  shown  in  Figure  23.  The  position  of 
moving  body  segment  1  relative  to  the  fixed  body  segment  2  is 
described  by  two  independent  coordinate  systems  shown  in 
Figure  23.  An  inertial  coordinate  system  (x,y)  with  unit 


Figure  23.  A  two-body  segmented  joint  is  illustrated  showing 
the  position  of  a  point,  Q,  attached  to  the  moving  coordinate 
system  (x ' ,y ' ) . 
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vectors  i  and  j  is  connected  to  the  fixed  body  segment,  while 
the  coordinate  system  (x',y')  with  unit  vectors  i*  and  j'  is 
attached  to  the  center  of  mass  of  the  moving  body  segment  1. 

The  motion  of  the  moving  (x',y')  system  relative  to  the 
fixed  (x,y)  system  may  be  characterized  by  three  quantities: 
the  translational  movement  of  the  origin  of  the  (x',y') 
system  in  the  x  and  y-directions ,  and  its  rotation,  a,  with 
respect  to  the  x  and  y-axis. 

Let  the  position  vector  of  the  origin  of  the  (x',y') 
system,  in  terms  of  the  fixed  system,  be  given  by: 


ro  = 


x  1  +  y„i 
o  1  oJ 


(5) 


Let  the  vector  Pq  be  the  position  vector  of  an  arbitrary  point 
Q,  on  the  moving  body  segment  in  the  base  (i',j').  Let  rn  be 

A.  *  \ 

the  position  vector  of  the  same  point  in  the  base  (i,j). 

That  is, 


pq  "  xQif  +  yqj’  ^ 

fQ  =  V  +  yQJ 

Referring  to  Figure  23,  for  vectors  p'  and  r,  the  following 
relationship  can  be  written: 

[rQ]  =  [rQ]  ♦  [T]  [ P q ]  (8) 

where  [T]  is  a  2x2  orthogonal  transformation  matrix.  The 
angular  orientation  of  the  (x',y')  system  with  respect  to  the 
(x,y)  system  is  specified  by  the  four  components  of  [T]  matrix 
Assuming  a  to  be  the  angle  between  the  positive  direction  of 
x-axis  and  the  positive  direction  of  x'-axis  (Figure  23),  then 
the  transformation  matrix  [T]  is  written  as: 


[T] 


cos  a  -sin  a 
sin  a  cos  a 


'(9) 
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Therefore,  the  position  in  the  (x,y)  system  of  any  given 
point  in  (x’,y')  system  can  be  determined  knowing  [T]  and  r  . 

CONTACT  CONDITIONS 

Assuming  rigid  body  contact  between  the  two  body  segments 
at  point  C  as  shown  in  Figure  24,  let  us  represent  the  contact 
surfaces  by  smooth  mathematical  functions  of  the  following 
form: 

Y  =  f2(x)  (10) 

y’  =  UD 

As  implied,  equations  (10)  and  (11)  represent  the  fixed  and 
moving  surfaces,  respectively. 


Figure  24.  A  two-body  segmented  joint  is  illustrated  showing 
contact  point,  C,  location  and  relevant  vectors. 


I 

The  position  vectors  of  the  contact  point,  C,  in  the 

A  A 

base  (i,j),  is  denoted  by: 


rc  =  xci  +  f1(xc)j  (12) 

and  the  corresponding  one  in  the  base  ( i ' » j  *  D •  is  given  by 

p'  -  x-i’  +  f2(x*)r  (13) 

Then  at  the  contact  point,  C,  the  following  relationship  must 
hold: 

[rc]  -  [rD]  +  [THp'1  (I4) 

This  is  a  part  of  the  geometric  compatibility  condition  for 
the  two  contacting  surfaces.  In  addition,  the  unit  normals 
to  the  surfaces  of  the  moving  and  fixed  body  segments  must  be 
colinear. 

A  A 

Let  n^  and  t^  be,  respectively,  the  unit  normal  and  unit 
tangential  vectors  to  the  fixed  surface,  y  =  f^x),  at  the 
contact  point,  C,  (Figure  25)  and  represented  by: 


nl  =  nlx1  +  nlyj 

h  =  ‘lx*  *  clyj 
or  for  t-^: 


t 


(15) 

(16) 

(17) 


By  definition,  a  unit  normal  vector  directed  toward  the  center 
of  curvature  is  given  by: 


»  dt, 

nl  =  R  35“ 


(18) 


where  R  is  the  radius  of  curvature  and  is  defined  as 
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SEGMENT  2 


Figure  25.  Unit  normal,  n^,  and  tangential,  t ^ ,  vectors  at 
the  contact  point,  C. 


[1  +  (df^dx)  2]3/2 

R  =  -  - = -  (19) 

|dZf1/dx2| 

The  distance,  ds,  along  the  surface,  y  =  f^Cx),  is  defined  as 
ds  =  y(dx)2  +  (dy)2 


or 


(20) 


Substituting  equation  (20)  in  equation  (17)  gives: 


tl 


1  +  (dfj/dx) 


Consequently,  equation  (18)  can  be  written  as: 


nl  = 


1  +  (dfj/dx) 


But  from  equation  (21) 


1  +  (dfj/dx) 


(df1/dx) (d2f1/dx2)  dr( 
[1  +  (dfj/dx)2]372  dx 


dtl  =  _ 1 _ 

[1  +  (dfj/dx) ] 1/2 


d2?  (df,/dx) (d2f,/dx2)  df 

— T  -  - 1 - S -  35T  (23) 

dxz  [1  +  (dfj/dxr]  ax 


Combining  equations  (22)  and  (23): 


nl  = 


[1  +  (dfj/dx) 2] 1/2 
|d2f1/dx2| 


(df-./dx)  (d^f-j/dx")  dr 
- i - -  -j-c  (24) 

[1  +  (dfj/dx)2]  Ux 


=  i  + 


Z5T  }  J 


■(*> 


Substituting  equations  (25)  and  (26)  into  equation  (24) 


[1  ♦  (df1/dx)2]1/2 
|d2f1/dx2| 


df,/dx  / 

j - j  ( 

[1  +  (dfj/dx)4-!  V 


i  +  -, 


‘c nr  J 
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or: 


«  (d  f^/dx  )  r  2  1  1/2 

II  =  — y-± - y-  1  +  (df,/dx)2  1/2 

1  |d2f1/dx2|  L  1  J 


dfx/dx 


[1  +  (dfj/dx)*] 


1  - 


(df^dx) 


[1  +  (df1/dx)A] 
Simplifying  the  above  equation  yields: 


nl  = 


(d2fx/dx2) 

<  \2  ’ 

-1/2  r 

/dM  ■  -1 

|d2f1/dx2| 

1  +  ( 

[' 

W  1  *  Jj 

where  x  is  the  x-coordinate  of  the  contact  point,  C,  in  the 
base  (i,j).  From  equations  (15)  and  (27)  it  can  be  shown  that 


2  2 
nlx  +  nly  "  1 


(28) 


which  is  an  expected  result  since  nlx  and  n^y  are  respectively 
the  x  and  y  components  of  unit  vector  n.  Similarly,  following 
the  same  procedure  as  was  outlined  for  n, ,  at  the  point  of 

A  1 

contact,  C,  the  unit  normal  vector,  n£ ,  to  the  surface, 
y'  =  f2(x')»  directed  toward  its  center  of  curvature,  can  be 
written  as: 


(29) 


where  x'Q  is  the  x 1 -coordinate  of  the  contact  point,  C,  in  the 
base  ( i ’  ,  j  * ) -  From  the  transformation  matrix  [T]  ,  given  in 
equation  (9),  the  coordinate  base  (i’,j')  is  related  to 

A  A 

coordinate  base  (i,j)  by: 


i'  *  cosa  i  +  sina  j 


(30) 
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j '  *  -sino  1  +  cosa  j 

Substituting  equations  (30)  and  (31)  into  equation  (29) , 
yields : 


(d2f2/dx  2) 

n?  =  — 1  + 
2  |d2f2/dx  2| 


/  df2  \  2 1  1/2  f  /  df2 

1  +  ( d7j  J  [’  (cosa  d7  + 


sina 


+  cosa  -  sina 


m 


The  colinearity  of  the  normals  at  the  point  of  contact,  C, 
requires  that: 

n^  x  n2  =  0  at  x=xc  »  x,=x' 

that  is: 

(ihi\/  |d2fi  //d2f2\/d2f2 

\  dx2  )'  ldx2  /  dfi  :  +  :  \  1  \  dx'2/  '  dx'2 

T  /  df,  \2  1/2  \  )  jT  /  df,  \  2  ‘ 

L1+  W)  J  <  L1  +  (57) 


[■< 


cosa  - j-  +  s 

dx 


ina^  i  +  ^  cosa  -  sina  — 2  ^  j  j 


dfj  df2  df^  df2 

sina  -r- —  - r  -  cosa  -r—  +  cosa  — r  +  sina  =  0 

dx  dx  dx  dx 

Finally,  the  contact  condition  takes  the  following  form: 


sina 


1  *  (^)x=xc  (£),.,c 


DYNAMIC  EQUATIONS  OF  THE  MOVING  BODY  SEGMENT 

The  total  load  acting  on  the  moving  body  segment  is 
divided  into  forces  exerted  by  nonlinear  springs  which  are 
simulating  the  ligament  and  capsule  forces,  contact  forces 
between  the  moving  and  fixed  body  segments  and  applied  exter¬ 
nal  forces  and  moments  (Figure  26).  In  the  following  sections 
each  of  these  forces  will  be  mathematically  formulated  and 
final  equations  of  motion  for  the  moving  body  segment  will  be 
presented. 

Ligament  Forces 

Representation  of  ligaments  and  capsules  as  nonlinear 
elastic  springs  along  with  their  governing  constitutive  equa¬ 
tions  were  discussed  previously.  The  force,  Fj ,  in  nonlinear 
spring,  j ,  is  a  function  of  its  length,  i . ,  that  is: 


Figure  26.  Forces  acting  on  the  moving  body  segment  of  a  two- 
body  segmented  joint. 


(37) 


Let  (p'O  be  the  position  vector  in  the  base  (i’,j')  of 
j  m  r 

the  insertion  point  of  the  ligament,  j,  in  the  moving  body 
segment.  The  position  vector  of  the  origin  point  of  the  same 
ligament,  j,  in  the  fixed  body  segment  is  denoted  by  (r2-)£ 

(S  A  J 

in  the  base  (i,j).  Here  the  subscripts  m  and  f  outside  the 
parenthesis  imply  "moving"  and  "fixed",  respectively.  The 
current  length  of  the  ligament  is  given  by: 


L. 


J 


(38) 


The  unit  vector,  X j ,  along  the  ligament,  j,  directed  from  the 
moving  to  the  fixed  body  segment  is: 


=  tj  [tf2j5 


Thus,  the  axial  force  in  the  ligament,  j,  in  its  vectorial 
form,  becomes: 


(39) 


Pj  -  Fjij  (40) 

where  Fj  is  given  by  equation  (3) . 

Contact  Forces 

Since  the  friction  force  between  the  moving  and  fixed 
body  segment  is  neglected,  the  contact  force  will  be  in  the 
direction  of  the  normal  to  the  surface  at  the  point  of 
contact.  The  contact  force,  N,  acting  on  the  moving  body 
segment  is  given  by: 

N  =  y  N  n^  at  x=xc  (41) 

where  N  is  the  magnitude  of  the  contact  force  and 
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(d2f1/dx2) 

|d2f1/dx2| 


at  x=x. 


Y  is  either  +1  or  -1  and  it  ensures  the  correct  direction  of 
the  contact  force  acting  on  the  moving  body  segment. 

Applied  External  Forces  and  Moments 

In  order  to  achieve  the  desired  joint  motion,  an  external 
force,  F  ,  and  a  moment,  M  ,  of  known  magnitude  and  direction 
are  applied  to  the  center  of  mass  of  the  moving  body  segment. 
The  force,  F  ,  and  moment,  M  ,  are  applied  in  the  base  (i,j) 
and  their  resultants  are  given  as: 

fe  *  <Vx  1  *  <Fe>y  3  (43) 

ft  -  M_  k  (44) 

e  e 

where  M£  is  the  magnitude  of  the  applied  moment  vector,  Mg . 
Equations  of  Motion 

The  dynamic  equations  of  motion  of  the  moving  body 
relative  to  the  fixed  body  segment  are  as  follows: 


x  =  Mxo 


(Fe)y  *  YN(nj)y  *  l  FjCijJy  -  My0 


M  *  (Tp  ' )  x  (yNii 


,)  ♦  l  (Tp!)  x  (F - Xj )  =  I,o 

1  j=l  J  33  Z 


where  p  is  the  number  of  ligaments  and  the  subscripts,  x  and 
y,  denote  the  components  of  the  related  quantities  in  the  x 
and  y-directions .  The  mass  of  the  moving  body  segment  is 
denoted  by  M  and  the  dots  denote  derivatives  with  respect  to 
time,  t.  The  mass  moment  of  inertia  of  the  moving  body  segment 
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about  the  z-axis  is  I  and  a  designates  its  angular 
acceleration.  The  problem  description  is  completed  by  assign¬ 
ing  the  initial  conditions,  which  are: 


Xn  "  Yn  = 


=  0 


(48) 


along  with  the  specified  values  for  xQ,  yQ  and  a  at  t  =  0. 

Three  nonlinear  second  order  differential  equations, 

(45) ,  (46)  and  (47)  ,  along  with  the  geometric  compatibility 
and  contact  conditions  of  equations  (14)  and  (36)  ,  provide 
the  necessary  relationships  in  order  to  determine  the 
following  unknowns: 

a)  xQ  and  yQ:  the  components  of  vector  rQ; 

b)  x  and  x ' :  the  x  and  x ' -coordinates  of  the 

^  ^  A  A  ^  A 

contact  point,  C,  in  the  base  (i,j)  and  (i',j')> 
respectively; 

c)  a:  the  orientation  angle  of  the  moving  (x',y') 
system  relative  to  the  fixed  (x,y)  system;  and 

d)  N:  the  magnitude  of  the  contact  force. 

The  numerical  procedure  employed  in  the  solution  of  the 
governing  equations  is  described  in  the  following  section. 

NUMERICAL  METHOD  OF  SOLUTION 

Newmark  Method  of  Differential  Approximation 

The  first  step  in  arriving  at  a  numerical  solution  of 
these  equations  is  the  replacement  of  the  time  derivatives 
with  a  temporal  operator;  in  the  present  work,  the  Newmark 
operators  (Bathe  and  Wilson  [1976])  are  chosen  for  this 
purpose.  For  instance,  xQ  is  expressed  in  the  following  form: 


_ 4 

Cat")1 


*r4t> 


4_  * t -  A t  "t-At 

At  X  Xo 


(49) 


•t-At 

o 


+ 


At  "t-At 
~ Z  xo 


(50) 


in  which  At  is  the  time  increment  and  the  superscripts  refer 

to  the  time  stations.  Similar  expressions  are  used  for  y  and 

o 
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a.  In  the  application  of  equations  (49)  and  (50),  the 
conditions  at  the  previous  time  station  (t-At)  are,  of  course, 
assumed  to  be  known. 

Numerical  Procedure 

After  the  time  derivatives  in  equations  (45)  ,  (46)  and 
(47)  are  replaced  with  the  temporal  operators  defined  in  the 
previous  section,  the  governing  equations  take  the  form  of  a 
set  of  nonlinear  algebraic  equations.  The  solution  of  these 
equations  is  accomplished  by  an  iteration  method.  In  this 
work,  the  Newton-Raphson  (Kao  [1974])  iteration  process  is 
used  for  the  solution.  To  linearize  the  resulting  set  of 
simultaneous  algebraic  equations,  we  assume: 


k  t  k-1  t  .  . 

+  Ax 

o  o  o 


(51) 


and  similar  expressions  for  the  other  variables  are  written. 
Here,  the  right  subscripts  denote  the  time  station  under 
consideration  and  the  left  subscripts  denote  the  iteration 
number.  At  each  iteration,  k,  the  values  of  the  variables  at 
the  previous  (k-1)  iteration  are  assumed  to  be  known.  The 
delta  quantities  denote  incremental  values.  Equation  (51)  and 
the  corresponding  ones  for  the  other  variables  are  substituted 
into  the  governing  nonlinear  algebraic  equations  and  the 
higher  order  terms  in  the  delta  quantities  are  dropped.  The 
set  of  n  simultaneous  algebraic  (now  linearized)  equations  can 
be  put  into  matrix  form: 

[K]  (A)  =  { D)  (52) 

where  [K]  is  an  nxn  coefficient  matrix,  {A}  is  a  vector  of 
incremental  quantities  and  {[)}  is  a  vector  of  known  values. 

The  iteration  process  at  a  fixed  time  station  continues 
until  the  delta  quantities  of  all  the  variables  become 
negligibly  small.  A  solution  is  accepted  and  the  iteration 
process  is  terminated  when  the  delta  quantities  become  less 
than  or  equal  to  a  prescribed  percentage  of  the  previous 


values  of  the  corresponding  variables.  The  converged  solutions 
of  each  variable  is  then  used  as  the  initial  value  for  the 
next  time  step  and  the  process  is  repeated  for  consecutive 
time  steps. 

The  only  problem  that  the  Newton -Raphson  process  may 
present  in  the  solution  of  dynamic  problems  is  due  to  the  fact 
that  the  period  of  the  forced  motion  of  the  system  may  turn 
out  to  be  quite  short.  In  this  case  it  becomes  necessary  to 
use  very  small  time  steps;  otherwise  a  significantly  large 
number  of  iterations  is  required  for  convergence.  This  matter 
will  be  further  dis:ussed  in  later  sections. 

Any  joint  can  be  modeled  with  this  general  formulation. 

The  knee  joint  will  be  modeled  accordingly  and  results  will 
be  presented  in  the  remaining  sections  of  this  report. 

MATHEMATICAL  DESCRIPTION  OP  A 
THREE-DIMENSIONAL  DYNAMIC  MODbL 

In  the  previous  section,  a  mathematical  description  for 
the  dynamic  motion  of  a  two-dimensional  articulating  joint 
was  presented.  Using  the  knowledge  and  the  insight  from  the 
discussion  along  with  the  mathematical  descriptions  of  the 
articular  surfaces,  ligaments  and  capsule  discussed  in  prior 
sections,  a  general  formulation  fo"  a  three-dimensional  mathe¬ 
matical  model  of  an  articulating  joint  will  be  presented  in 
this  sect  ion . 

Once  again  for  the  purpose  of  studying  the  joint  section, 
one  segment  is  assumed  to  be  fixed  while  the  other  segment  is 
executing  a  relative  motion.  Coordinate  systems  (x,y,z)  and 
(x',y',z')  are  attached  to  the  fixed  and  moving  body  segments, 
respectively,  and  their  relative  position  is  discussed  below. 
This  relative  position  is  determined  by  six  variables:  three 
components  of  a  vector  specifying  the  origin  of  the  (x',y',r') 
system,  and  three  rotations  which  determine  the  orientation  of 
the  (x',y',z')  system.  Contact  conditions  and  dynamic  equation 
of  motion  are  then  presented. 


CHARACTERIZATION  OF  RELATIVE  POSITIONS 

A  joint  connects  two  body  members  which  are  here 
designated  as  segment  1  and  segment  2.  The  position  of  the 
moving  body  segment  1  relative  to  fixed  body  segment  2  is 
described  by  two  coordinate  systems  as  shown  in  Figure  27. 
Inertial  coordinate  system  (x,y,z)  with  unit  vectors  i,j  and 

A 

k  is  connected  to  the  fixed  body  segment  and  coordinate  system 

A  A  A 

(x',y',z')  with  unit  vectors  i’,j'  and  k'  is  attached  to  the 
center  of  mass  of  the  moving  body  segment.  The  (x',y',z') 
coordinate  system  is  also  taken  to  be  the  principal  axis 
system  of  the  moving  body  segment. 

The  motion  of  the  moving  (x',y',z')  system  relative  to 
the  fixed  (x,y,z)  system  may  be  characterized  by  six 
quantities:  the  translational  movement  of  the  origin  of  the 


Figure  27.  A  two-body  segmented  joint  is  illustrated  in 
three  dimensions,  showing  the  position  of  a  point,  Q,  attached 
to  the  moving  coordinate  system  (x',y',z'). 


(x',y',z')  system  in  the  x,  y  and  z -directions ,  and  6,  <J>  and  4* 
rotations  with  respect  to  the  x,  y,  and  z-axes. 

Let  the  position  vector  of  the  origin  of  the  (x',y',z') 
system  in  the  fixed  system  be  given  by  (Figure  27) : 


r  =  xi  + 
o  o 


+  zok 


(53) 


Let  the  vector,  p ’ ,  be  the  position  vector  of  an  arbitrary 

^  AAA 

point,  Q,  on  the  moving  body  segment  in  the  base  (i',j’,k'). 
Let  rq  be  the  position  vector  of  the  same  point  in  the  base 
(i, j ,k) .  That  is. 


(54) 


AAA 

Tq  =  xQi  +  ypj  +  Zgk  (55) 

Referring  to  Figure  27,  vectors  Pq  and  fp  have  the  relation¬ 
ship  : 

{rq}  =  {ro}  +  fTftpq}  (56) 

where  [T]  is  a  3x3  orthogonal  transformation  matrix.  The 
angular  orientation  of  the  (x',y',z')  system  with  respect  to 
the  (x,y,z)  system  is  specified  by  the  nine  components  of  the 
[T]  matrix  and  can  be  written  as  a  function  of  the  three 
variables,  e,  <p  and  \p : 


T  =  T(e,<t>,40  (57) 

There  are  several  systems  of  variables  such  as  e ,  <fi  and 
ip  which  can  be  used  to  specify  T.  In  this  study  the  Euler 
angles  will  be  utilized. 

The  orientation  of  the  moving  coordinate  system 
(i',j’,k')  is  obtained  from  the  fixed  coordinate  system 
(i,j,k)  by  applying  successive  rotation  angles,  <j> ,  6  and  4> 
(Figure  28).  First  the  (i,j,k)  system  is  rotated  through  an 
angle  $  about  the  z-axis  (Figure  28a)  ,  which  results  in  the 
intermediary  system  (i^,jj,kj),  where: 
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Figure  28.  Successive  rotations  of  <j> ,  e  and  4 /,  of  the  (x,y,z) 
coordinate  system. 


h 

=  cosi(i  i  +  s in4>  j 

(58a) 

h 

=  -  sin<j>  i  +  cos4>  j 

(58b) 

*1 

=  k 

(58c) 

The 

(Figure 

where : 

second  rotation  through  an  angle  e  about  the 
28b),  produces  the  intermediary  system  C i 2 » j  2 

i^ -axis 

* 

h 

=  *1 

(59a) 

h 

=  cos6  j  j  +  sine 

(59b) 

^2 

=  -  sine  jj  +  cose  k^ 

(59c) 
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The  third  rotation  through  an  angle  ip  about  the  k2*axis 
(Figure  28c),  gives  the  final  orientation  of  the  moving 

AAA  A  A  A 

(i',j',k')  system  relative  to  the  fixed  (i,j,k)  system,  where 


i  ’  =  1  ^  =  cosiji  i 2  +  sinij/  j  2 
j '  =  j  2  =  -  sinijj  i2  +  costj;  j  2 


k'  -  k3  k2 


Substituting  equations  (58)  and  (59)  into  equation  (60),  the 
final  orientation  of  the  (i',j',k')  system  relative  to  the 
(i,j,k)  system  may  be  written  as: 

i'  =  [  (cos4>  cos<}>  -  sini|/  cose  sintj>)i  +  (cosip  sin<}> 

+  sirnji  cose  cos<f>)j  +  (simj)  sine)k] 

j'  =  [(-simji  coscj)  -  cosi|>  cose  sin<J>)i  +  (-  sinijj  sin4> 

+  cosij>  cos <j>  cose)  j  +  (costjj  sine)k] 

k'  =  [(sine  sin<())i  -  (sine  cos$)j  +  cose  k] 
or,  in  matrix  notation: 


r 

k' 
where 
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CONTACT  CONDITIONS 

Assuming  rigid  body  contacts  between  the  two  body 
segments  at  points  C^  (i=l,2)  as  shown  in  Figure  27,  let  us 
represent  the  contact  surfaces  by  smooth  mathematical  func¬ 
tions  of  the  form: 

y  =  f(x,z)  (64) 

y'  =  g(x’  ,z')  (65) 

As  implied,  equations  (64)  and  (65)  represent  the  fixed  and 
the  moving  surfaces,  respectively. 

The  position  vectors  of  the  contact  points  C-  (i=l,2)  in 

AAA  J- 

the  base  (i,j,k)  is  denoted  by 

rc  =  x  i  +  f(x  ,  z  )j  +  z  k  (66) 

i  ci  ct  Ci 

and  the  corresponding  ones  in  the  base  ( i ' , j  r , k * )  are  given 
by: 


p 


i 


+ 


+ 


(67) 


Then,  at  each  contact  point  C^  the  following  relationship  must 
hold: 


{rCi}  =  {ro}  +  lTJT{pc.}  (68) 

This  is  a  part  of  the  geometric  compatibility  condition  for 
the  two  contacting  surfaces.  Furthermore,  the  unit  normals 
to  the  surfaces  of  the  moving  and  fixed  body  segments  at  the 
points  of  contacts  must  be  colinear. 

Let  hc  ( i  =  1 , 2)  be  the  unit  normals  to  the  fixed  surface, 
y  *  f(x,z),1at  the  contact  points,  ( i* 1,2),  then: 


where  rc^  is  given  in  equation  (66)  and  the  components  of  the 
matrix  [G]  are  determined  by: 
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Since  (3zc^/3xc^)  =  0  and  (3xc^/3zc^)  =  0,  then  the  components 
of  matrix  [GUol  reduce  to: 


From  equations  (72),  the  DET[G]  can  be  written  as: 


DET [G]  =  1  + 


and  therefore,  the  unit  outward  normals  expressed  in  equation 
(69)  will  have  the  following  form: 


i  -  i 


where  the  parameter,  y,  is  chosen  such  that  n  represents  the 

c  i 

outward  normal. 

Similarly,  following  the  same  procedure  as  outlined  above 
n£^  (i=l,2),  the  unit  outward  normal  to  the  moving  surface, 
y'  =  g(x',z'),  at  contact  points,  C-  (i=l,2),  and  expressed  in 

AAA  ^ 

( i *  , j '  , k ' )  system,  can  be  written  as: 


1  ♦  f  ^T- 

I  5  V  ' 


■3g 


^r'-] 


a  7  ' 


where  parameter,  8,  is  chosen  such  that  n^  represents  the 
outward  normal. 

Colinarity  of  unit  normals  at  each  contact  point 
(i=l,2),  requires  that: 

(n  }  =  -[T]T{n'  } 

Li  i 

Note  that  colinearity  condition  can  also  be  satisfied  by 

T 

requiring  that  the  cross  product  (nc^  x  T  n^)  be  zero. 

DYNAMIC  EQUATIONS  OF  THE  MOVING  BODY  SEGMENT 

The  total  load  acting  on  the  moving  body  segment  is 
divided  into  forces  exerted  by  nonlinear  springs  which 
simulate  the  ligaments  and  capsule  forces,  contact  forces 


1 


between  the  moving  and  the  fixed  body  segments  and  applied 
external  forces  and  moments  (Figure  29).  Ligament  and 
capsule  forces  were  discussed  and  formulated  previously.  The 
coefficient  of  friction  between  the  articulating  surfaces, 
owing  to  the  presence  of  the  synovial  fluid  in  the  joints,  is 
known  to  be  very  low  (Radin  and  Paul  [1972]).  Accordingly, 
the  friction  force  between  the  two  body  segments  will  be 
neglected.  Therefore,  the  contact  forces,  N^,  acting  on  the 
moving  body  segment  are  given  by: 


=  I  Ni  I  f(nc.)xi  +  (nc.M  *  (nc.)Zk] 


(77) 


where  |N^|'s  are  the  unknown  magnitudes  of  the  contact  forces 
and  (nc^)x,  (nc^)y  and  (nc^)z  are  the  components  of  the  unit 
normal,  nc^,  in  the  x,  y  and  z-directions .  respectively. 


Figure  29.  Forces  acting  on  the  moving  body  segment  of  a  two- 
body  segmented  joint  in  three  dimensions. 
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The  desired  dynamic  motion  is  achieved  by  applying  an 
external  force  and  moment  to  the  center  of  mass  of  the  moving 
body  segment  whose  resultants  are  given  as: 

Fe  -  (Fe)xl  +  (Fe)yj  ♦  (Fe)zk  (78) 

«e  -  (Me)xi  *  (Me)yj  ♦  CMe) 2k  (79) 

Equations  of  motion  will  be  presented  next. 

Equations  of  Motion 

The  equations  governing  the  forced  motion  of  the  moving 
body  segment  are: 

q  P 


(Fe>x  *  Jj  lNil<nci)x  *  Fj  (xj '  x  ■  Mxo 

(Fe’y  *  jj  *  jj  Wy  '  < 

(Fe>z  *  Jj  |NtIUCi)z  *  FjUj)z  *  M'z0 


(80a) 

(80b) 

(80c) 


x'  x' 

X 
'  3 

X 

X 

p— < 

II 

+  z ' z  ' 

ly'y' 

)uiy,a)z, 

(81a) 

y’y' 

=  Iy'y'“y' 

+  ^ 1 X  *  X  ' 

1  z '  z  ' 

)coz,u>  - 

(81b) 

z '  z ' 

*  z ’ z ’ “z ' 

+  ^y'y' 

)wx* “y* 

(81c) 

where  p  and  q  represent  the  number  of  ligaments  and  the 
contact  point,  respectively.  Ix,x,f  I  ,  ,  and  Iz,2,  are  the 
principal  moments  of  inertia  of  the  moving  body  segment  about 
its  centroidal  principal  axis  system  (x',y',z'),  and  w  , ,  u  , 
and  u>z,,  are  the  components  of  the  angular  velocity  vector 
which  are  given  below  in  terms  of  the  Euler  angles: 


=  0  cosi|i  +  |  sine  sinifi 


(82a) 


u)y,  =  -0  sintj/  +  if  sine  cosij* 
a)  ,  =  |  cose  +  4> 

4 

The  angular  acceleration  components,  u»x ,  u>  ,  and  uz  are 
directly  obtained  from  equation  (82): 

<Lx,  =  e  costjj  -  ( 0  sin;ji  -  $  cosi|»  sine) 

+  <j>  sine  sin^  +  <Je  cose  sin^ 

o)y,  =  -  e  sim()  -  4> ( ©  cosi{/  +  sin^  sine) 

+  <}>  sine  cost}/  +  <Je  cose  cos\j> 


(82b) 

(82c) 

(83a) 

(83b) 


•  '*  •  •  ** 

ojz,  =  <j>  cose  -  <j>e  sine  +  ip  (83c) 

Note  that  the  moments  components  shown  on  the  left-hand  side 
of  equation  81  are  obtained  from: 

q  p 

M  =  Me  +  l  [T]T(p'  )  x  (|N.|n  )  +  l  [T]T(p '-)  x  ( F  -  X  . )  (84) 

i=l  i  1  i  j=X  J  J  J 

where  Mg  is  applied  external  moment,  and  p  and  q  represent 
the  number  of  ligaments  and  contact  points,  respectively. 

The  equations  (80)  and  (81)  form  a  set  of  six  nonlinear 
second  order  differential  equations  which,  together  with  the 
contact  conditions  (68)  and  (76)  ,  form  a  set  of  16  nonlinear 
equations  (assuming  two  contact  points,  i.e.  i*l,2)  with  16 
unknowns : 

a)  d ,  <p  and  ip ,  which  determine  the  components  of 
transformation  matrix  [T] ; 

b)  x  ,  yQ,  and  zQ:  the  components  of  position  vector 


c)  xc^,  zc^,  x^  and  z^  (i  =  l,2):  the  coordinates  of 
contact  points; 

d)  | N ( i= 1,2):  the  magnitudes  of  the  contact  forces. 


Numerical  procedure  outlined  previously  can  be  utilized 
for  the  solution  of  the  three-dimensional  joint  model 
equations  presented  in  this  section.  However,  becuase  of  the 
extreme  complexity  of  these  equations,  in  this  report  we  will 
present  in  some  detail  only  the  numerical  solution  of  a  two- 
dimensional  joint  model  applied  to  the  human  knee  joint. 

TWO-DIMENSIONAL  DYNAMIC  MODEL 
OF  THE  KNEE  JOINT 

Thus  far,  mathematical  descriptions  of  the  articular 
surfaces,  the  ligaments  and  capsule  have  been  developed  as 
well  as  general  formulations  for  two-  and  three-dimensional 
dynamic  model  of  an  articulating  joint.  These  formulations 
can  now  be  applied  to  a  mathematical  description  for  the 
dynamic  motion  of  the  knee  joint.  The  most  general  and 
realistic  model  of  the  knee  should  be  three-dimensional. 
However,  a  simpler  two-dimensional  model  can  be  helpful  and 
rewarding  in  understanding  the  essentials  of  the  problem  and 
serve  as  the  groundwork  for  the  sound  development  of  a  three- 
dimensional  model.  Before  we  present  our  two-dimensional 
dynamic  model  of  the  knee  joint  we  will  briefly  discuss  the 
previously  developed  models  of  the  knee. 

PREVIOUSLY  DEVELOPED  KNEE  JOINT  MODELS 

The  models  of  the  knee  joint  can  be  subdivided  into  pure 
kinematical  models  and  models  describing  the  force  action  in 
the  joint.  Kinematical  models  try  to  describe  the  motions 
between  femur  and  tibia  without  considering  the  forces  and 
moments  in  the  joint.  A  model  of  this  type,  developed  by 
Strasser  [1917]  is  a  four-bar  mechanism  (Figure  30).  Two 
bars  represent  the  cruciates,  while  the  other  bars  represent 
femur  and  tibia,  respectively.  Menschik  [ 1974a  ,  1974b] 
extended  this  planar  model  by  two  curves  representing  the 
tibial  and  femoral  articular  surfaces  and  also  studied 
location  of  the  insertion  areas  of  the  collateral  ligaments 
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Figure  30.  Four-bar  mechanism  according  to  Strasser  [1917] 
is  shown.  Components  are  identified  as  1  =  tibia,  2  = 
anterior  cruciate  ligament,  3  =  femur  and  4  =  posterior 
cruciate  ligament. 


in  this  model.  Similarly,  the  model  of  Huson  [1974,1976]  is 
based  on  the  idea  of  a  four-bar  mechanism,  but  it  is  also 
able  to  simulate  internal -external  free  range  of  motion  by 
means  of  a  certain  inclination  of  the  plane  representing  the 
lateral  tibial  articular  surface. 

Several  planar  mechanisms  simulating  the  motion  of  the 
human  knee  joint  in  the  sagittal  plane  have  been  proposed  by 
Freudenstein  and  Woo  [1969].  The  aim  of  this  study  was  to 
serve  as  a  guideline  in  the  kinematical  design  of  joint 
prostheses.  Investigations  aimed  at  accurate  in-vivo  measure 
ments  of  three-dimensional  relative  motions  in  a  human  knee 
joint  so  far  have  not  been  fully  described  in  the  literature. 
A  study  by  Levens ,  Inman  and  Blosser  [1948]  was  limited  to 
relative  rotations  in  a  transverse  plane.  In  this  study, 
stainless  steel  pins  of  2.5  mm  diameter  were  drilled  firmly 
into  the  femur  and  the  tibia,  sterility  precautions  and  local 
anesthesia  being  used.  No  further  attempts  using  this  method 
are  described  in  the  literature. 

The  mechanical  analysis  of  the  human  knee  joint  has  in 
the  past  been  carried  out  mostly  with  human  knee  joint 


specimens  which  can  be  considered  as  the  best  available 
representation  of  the  living  human  joint.  The  shortcomings 
of  the  knee  specimens  are  the  lack  of  muscle  forces  and  the 
difference  in  material  properties  from  those  of  the  living 
joint.  Moreover,  the  availability  of  the  knee  specimens  for 
research  purposes  is  limited.  Living  and  dead  animals  have 
often  been  used  as  surrogates  for  human  material  in  the  study 
of  special  problems.  For  example,  the  mechanical  behavior  of 
the  menisci  in  canines  and  pigs  have  been  studied  by  Krause 
et  al.  [1976]  and  Jaspers  et  al.  [1978],  respectively. 

In  general,  kinematical  models  offer  valuable 
possibilities  of  gaining  a  better  insight  into  some  aspects 
of  joint  behavior.  However,  their  application  is  restricted 
to  phenomena  in  which  force  actions  are  of  no  interest. 

A  number  of  knee  models,  described  in  literature,  were 
not  concerned  with  relative  joint  motions  but  were  developed 
to  determine  the  forces  in  ligaments,  muscles  or  between 
articular  surfaces.  In  these  models  the  joint  structures  are 
simplified  in  such  a  way  that  a  so-called  statically  determi¬ 
nant  system  is  achieved.  Consequently,  the  constitutive 
equations  of  the  ligaments  are  not  necessary  and  the  condi¬ 
tions  of  equilibrium  are  sufficient  to  determine  the  relevant 
forces.  Two-dimensional  models  of  this  type  are  reported  by 
Kettelkamp  and  Chao  [1972],  Smidt  [1973],  Perry,  Antonelli  and 
Ford  [1975],  Seedhom  and  Terayama  [1976],  Rens  and  Huiskes 
[1976],  while  Morrison  [1967,1970]  presented  a  three- 
dimensional  model. 

Although  the  studies  on  the  biomechanics  of  the  knee 
joint  have  a  long  history,  those  studies  which  are  essentially 
a  mathematical  modeling  of  this  largest  and,  apparently,  most 
complicated  joint  in  the  body  are  few.  Crowninshield,  Pope 
and  Johnson  [1976]  felt  justified  in  stating  that  there  were, 
at  that  time,  no  analytical  models  of  the  knee  available 
which  permit  the  prediction  of  the  response  of  the  joint  to 
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either  external  forces  or  displacements.  In  their  analytical 
model,  on  the  other  hand,  only  the  quasi-static  response  of 
the  joint  was  studied  and  the  overall  joint  stiffness  was 
obtained  as  a  function  of  the  flexion  angle.  The  method  used 
by  Crowninshield,  Pope  and  Johnson  [1976]  is  the  so-called 
inverse  method  in  which  the  ligament  forces  caused  by  a 
specified  set  of  translations  and  rotations  in  specified 
directions  are  determined  by  comparing  the  geometries  of  the 
initial  and  displaced  configurations  of  the  knee  joint.  In 
this  method  the  externally  applied  displacements  do  not  need 
to  be  continuous;  however,  the  discrete  values  used  are  to  be 
realistic.  In  Crowninshield,  Pope  and  Johnson  [1976]  these 
values  were  based  on  experimental  data  available  in  the 
literature.  The  purpose  of  the  work  was  to  obtain  the  stiff¬ 
ness  of  the  joint  which,  in  turn,  required  the  calculation  of 
the  ligament  lengths  at  various  knee  configurations.  Thus, 
it  was  not  necessary  to  consider  the  contribution  of  the 
curved  joint  surfaces  to  the  overall  mechanical  behavior  of 
the  knee.  Moreover,  the  external  forces  required  for 
equilibrium  were  not  determined  either. 

Improvements  to  the  quasi-static  model  discussed  above 
have  been  provided  recently  by  Wismans  et  al.  [1980],  in 
which  a  three-dimensional  analytical  model  of  the  femoro- 
tibial  joint  is  presented.  The  study  considers  not  only  the 
geometry  but  also  the  static  equilibrium  of  the  system.  The 
three-dimensional  curved  geometry  of  the  joint  surfaces  are 
included  in  the  model.  Ligaments  are  modeled  as  non-linear 
elastic  springs.  The  solution  method  employed  by  Wismans 
et  al.  [1980]  is  also  a  quasi-static,  inverse  method.  The 
flexion-extension  motion  is  simulated  by  prescribing  several 
flexion-extension  angles.  The  dependent  variables  of  the 
problem,  including  the  ligament  forces,  are  determined  from 
the  equilibrium  equations  and  the  geometric  compatibility 
conditions.  However,  for  nonlinear  problems  of  this  kind  it 
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is  known  that  there  can  exist  more  than  one  equilibrium 
configuration  for  a  given  flexion-extension  angle  unless  an 
external  force  is  also  specified.  Accordingly,  in  the 
inverse  method  utilized  by  Wismans  et  al.  [1980]  it  is 
necessary  to  specify  the  external  force  required  for  the 
preferred  equilibrium  configuration.  Such  an  approach  is 
applicable  only  in  quasi-static  analysis.  In  dynamic  analysis, 
the  equilibrium  configuration  preferred  by  the  system  is  the 
unknown  and  the  mathematical  analysis  itself  is  to  provide 
that  equilibrium  configuration. 

Biomechanics  of  the  knee  joint  has  also  been  investigated 
by  Andriacchi  et  al.  [1977,1978].  They  reported  a  statically 
indeterminate  model  for  the  analysis  of  motion  and  forces  in 
the  knee  joint.  Like  Crowninshield ,  Pope  and  Johnson  [1976], 
they  represent  ligaments  and  capsule  by  a  number  of  springs, 
while  the  articular  surfaces  and  menisci  are  also  represented 
in  the  model.  Numerical  predictions  are  consistent  with 
experimental  observations.  The  models  of  Andriacchi  et  al. 
[1977,1978]  and  Wismans  et  al.  [1980],  are  essentially  the 
same  in  the  sense  that  both  deal  with  the  quasi-static 
response  of  the  knee  joint.  Detailed  discussions  of  various 
anatomical  and  functional  aspects  of  the  human  knee  joint  can 
also  be  found  in  Gray  [1973],  Engin  and  Korde  [1974], 
Blacharski,  Somerset  and  Murray  [1975],  Jacobsen  [1976],  Pope, 
et  al.  [1976],  and  Engin  [1978]. 

As  seen  from  the  preceding  discussion,  mathematical 
modeling  of  the  knee  joint  has  not  yet  reached  a  definitive 
stage  of  development.  It  is  interesting  to  note  that  a 
biodynamic  model  of  the  knee  joint  is,  to  the  best  of  the 
authors'  knowledge,  not  yet  available  in  the  literature.  It 
is  more  appropriate  to  study  via  a  dynamic  model  the  response 
of  the  joint  to  dynamically  applied  loads.  The  artificial 
restrictions  of  the  quasi-static  inverse  method,  such  as  the 
necessity  to  specify  the  preferred  configuration,  can  be 


eliminated  if  the  dynamics  of  the  problem  are  incorporated 
into  the  model. 

TWO-DIMENSIONAL  GEOMETRY  OF  THE  KNEE  JOINT 
Coordinate  Systems 

Previous  sections  have  discussed  the  relative  positions 
and  coordinate  systems  for  a  general,  two-body  segmented 
joint.  Uniquely  specifying  the  locations  of  these  coordinate 
systems  is  virtually  impossible  because  of  the  lack  of  well- 
defined  anatomical  landmarks  in  the  human  body.  However,  the 
relative  motions  in  a  joint  are  not  affected  by  the  choice  of 
coordinate  systems  and  only  some  parameters  describing  the 
relative  displacements  have  different  values  for  different 
coordinate  systems. 

The  position  of  the  tibia,  defined  as  the  moving  segment, 
relative  to  the  femur,  defined  as  the  fixed  segment,  is  shown 
in  Figure  31.  The  origin  of  the  moving  coordinate  system 
(x 1 ,y ’ )  coincides  with  the  center  of  mass  of  the  tibia,  with 
the  y'-axis  being  directed  along  the  longitudinal  axis  of  the 
tibia.  The  inertial  coordinate  system  (x,y)  is  attached  to 
the  fixed  femur  with  the  x-axis  directed  along  the  posterior- 
anterior  direction  and  the  y-axis  coinciding  with  the  femoral 
longitudinal  axis.  The  locations  of  the  origins  of  these 
coordinate  systems,  shown  in  Figure  31,  are  4.01  cm  and 
21.34  cm  from  the  intersection  points  of  the  y  and  y'  axes 
with  the  articulating  surfaces.  The  rotation  of  the  moving 
(x',yf)  coordinate  system  with  respect  to  the  fixed  (x,y) 
system  is  denoted  by  a. 

Mathematical  Descriptions  of  the  Articulating  Surfaces 

The  two  dimensional  profiles  in  the  plane  of  motion  of 
the  femoral  and  tibial  articulating  surfaces  are  obtained 
from  X-ray  of  a  human  knee  joint.  The  coordinates  (x]<»yjc) 
of  a  number  of  points  on  these  profiles  are  measured  using  a 
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Figure  31.  Coordinate  systems  locations  and  relative  positions 
of  the  tibia  and  femur  are  shown  for  the  two-dimensional 
dynamic  model  of  the  knee  joint. 


two-dimensional  sonic  digitizing  technique  described  previously 
(Figure  32) . 

A  polynomial  equation  of  degree  n(n>l)  is  used  as  an 
approximate  mathematical  representation  of  the  profile  under 
consideration.  This  polynomial  has  the  form 

n 

y (x)  *  l  Aix1  (85) 

l 

where  A^'s  are  the  polynomial  coefficients.  These  coefficients 
are  determined  by  the  use  of  the  subroutine  program  CHEPLS 
(Appendix  A).  This  subroutine  is  capable  of  determining,  by 
means  of  statistical  tests,  where  the  set  of  given  data  points 
are  linear  or  nonlinear.  If  the  data  are  nonlinear  at  the 
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Stylus  Menu 


Figure  32.  Two-dimensional  microphone/sensor  with  menu 
capability,  and  sonic  emitter/stylus  for  the  Graf/Pen  sonic 
digitizer  is  shown  with  x-ray  of  the  knee  joint  in  place  for 
obtaining  the  profiles  of  the  articulating  surfaces. 


95%  confidence  level,  then  the  routine  finds  the  lowest 
degree  polynomial  which  adequately  represents  the  data.  The 
calculation  procedure  is  to  compare  the  standard  deviation  of 
the  model  of  degree  n  with  that  of  (n+l)**1  degree  model.  If 
there  is  no  significant  difference  at  the  95%  confidence 
level,  then  the  polynomial  of  degree  n  is  accepted  as  the 
lowest  degree  polynomial  which  adequately  represents  the  data. 
In  the  present  study,  this  method  yields  the  following 
equations  for  the  femoral  and  tibial  profiles,  respectively: 

f,(x)  =  0.04014  -  0.247621  x  -  6.889185  x2  -  270.4456  x3 

4  (86) 

-8589.942  x4 
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f2(x')  =  0.213373  -  0.0456051  x'  +  1.073446  x' 


(87 


with  the  aid  of  the  Versatec  plotter,  the  tibial  and  femoral 
articulating  surfaces  predicted  by  equations  (86)  and  (87) 
are  shown  in  Figure  33. 

DESCRIPTION  OF  THE  LIGAMENT  MODEL 

Selection  of  the  Springs  and  Corresponding  Parameters 

Only  four  major  ligaments  of  the  knee  joint  will  be 
considered  in  the  present  work  although  consideration  of  any 
other  ligament  presents  no  difficulty.  These  ligaments  are 
the  lateral  collateral  (LC) ,  the  medial  collateral  (MC) ,  the 
anterior  cruciate  (AC) ,  and  the  posterior  cruciate  (PC) .  Th 
ligaments  are  modeled  as  nonlinear  elastic  springs  having  a 


Figure  33.  Versatec  plot  of  the  two-dimensional  representa¬ 
tions  of  the  tibial  and  femoral  articulating  surfaces. 


constitutive  relation  given  by  equation  (3).  The  stiffness 
values  for  these  ligaments  are  summarized  in  Table  8.  Initial 
strains  in  the  ligaments  are  taken  as  zero  since,  at  present, 
there  is  no  accurate  data  available  on  strains  as  a  function 
or  flexion  angle.  Zero  strain  condition  for  the  ligaments 
can  be  partially  justified  if  an  appropriate  starting  flexion 
angle  under  no  external  load  is  chosen. 

Insertions  and  Origins 

The  coordinates  of  the  insertion  points  of  the  ligaments 
in  both  the  tibia  and  the  femur  are  determined  from  the 
available  information  in  the  literature  (Wang,  Walker  and 
Wolf,  [1973],  Crowninshield,  Pope  and  Johnson  [1976])  and  close 
anatomical  study  of  the  knee  joint.  These  coordinates  on  the 
tibia  are  denoted  by  xj  and  yt ,  and  those  on  the  femur  are 
denoted  by  Xj  and  y  j .  The  values  used  in  the  present  work 
are  summarized  in  Table  10.  Obviously,  these  values  are 
determined  with  respect  to  coordinate  systems  shown  in  Figure 
31,  and  for  the  specific  specimen  used  in  our  study. 


Table  10 


COORDINATE  VALUES  FOR  THE  INSERTIONS  AND  ORIGINS 
OF  THE  KNEE  JOINT  LIGAMENTS,  IN  METERS. 


LIGAMENT 

TIBIA 

FEMUR 

x! 

1 

y- 

X3 

Medial  Collateral 
(MC) 

0.008 

0.163 

-0.023 

0.014 

Lateral  Collateral 
(LC) 

0.025 

0.178 

-0.025 

0.019 

Anterior  Cruciate 
(AC) 

-0.005 

0.213- 

-0.023 

0.019 

Posterior  Cruciate 
(PC) 

0.025 

0.208 

-0.032 

0.024 

MATHEMATICAL  DESCRIPTION  OF  THE  DYNAMIC 
MOTION  OF  THE  KNEE  JOINT 

The  total  external  forces  acting  on  the  tibia  are  shown 
in  Figure  34.  These  forces  are:  ligament  forces,  Fj 
(j=l,--4);  normal  contact  force  N;  and  the  applied  external 
force,  F  ,  and  moment,  M  .  Frictional  forces  between  the 
femoral  and  tibial  surfaces  will  be  neglected  since  the 
coefficient  of  friction  between  the  articulating  surfaces, 
owing  to  the  presence  of  the  synovial  fluid,  is  known  to  be 
very  low  (Radin  and  Paul  [1972];  Dowson  [1976]).  Therefore, 
the  contact  force,  N,  will  be  in  the  direction  of  the  normal 
to  the  surface  at  the  point  of  contact. 

The  equations  governing  the  forced  motion  of  the  tibia 

are : 


Figure  34.  Forces  acting  on  the  moving  tibia  are  shown  for 
the  two-dimensional  model  of  the  knee  joint. 


CFe)x  *  vN(n1)x  * 
(Fe)y  *  YN(ni)y  -  \ 
Me  +  (Tp^.)  x  (yNn^  + 


Wx  =  Mxo 


Fj(Xj)y  =My0 


I  (Tpl)  x  (F.X  ) 
j=l  3  33 


V 


(88) 

(89) 

(90) 


where  all  the  parameters  in  the  above  governing  equations  have 
been  defined  in  previous  discussion  on  the  two-dimensional 
general  formulation  for  dynamic  motion. 

NUMERICAL  ANALYSIS  AND  DERIVATIONS 

The  differential  variables,  xQ,  yQ,  and  a  are  substituted 
by  their  corresponding  Newmark  approximations  given  in  equa¬ 
tions  (49)  and  (50).  Subsequently,  these  simultaneous  non¬ 
linear  algebraic  equations  along  with  geometric  compatibility 
and  contact  conditions  of  equations  (8)  and  (36)  are 
linearized  by  applying  the  Newton-Raphson  iteration  process. 
After  considerable  amount  of  mathematical  operations  and 
eliminations  of  the  second  order  variational  terms  (see 
Appendix  B  for  complete  derivations) ,  the  following  final  form 
of  the  linearized  governing  equations  of  motion  are  obtained: 


(yI11x)6N  +  (ynK)  v  Z'-x 

iX  w  nlx  At  xo 


(-%«, 


x  -  CVx 


-yNkn^x  +  M  ^  — — 9  [xk  -  x 


*  1  Fi 

j=l  3 

(4k  .  t-At,  .  4,-t-At  .  ”t-At j 

I  At^  0  0  Ato  O  I 


*  (yNk)6nly  - 


£  Fjy  •  CFe’y 


vNknky  ♦  M  j^jtyJ  -  y0 


t-At,  .  _4  •  t-At  .  "  t-At  ) 
1  At  yO  ' O  j 
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y  f.  6 


£  V*j  '  *  yN  ^  ,ii  “J*>J 

♦  [  j  fjx  ♦  vNkn?  ]J  M(xk  -  xk)nk  -  (yk  -  yk)nkx]«N 
j  =  1  J  '  o 

k_k 


-  ^NMx>5yc 


4  -t-At 
At  ° 


”\x 

+ 

41, 

«o 

4  rg 

4-> 

< 

1 

a 

- 

t-At  )  _ 

M 

.rkr„k 
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V  rx^  -  xk)F.  +  T  (yk  -  yk)F- 
j  =  lL  J  ojrjy  J  Y°  J* 

k .  k  I  ,  T  (  4  ,  k  t-At. 

y05nix>  +  I^^2ta  ) 


i 

Similarly,  the  linearized  geometric  constraint  equations  are: 
+  (x'k  sinak  +  y'k  cosak)6  -  cosak6  , 

C  C  Ot  a 


6  -  6 
X^_  x„ 

c  o 


k  k 

+  sina  6  ,  =  x£  -  X 
yc 


I 

'c 


k  .  v.k  k  . k  k 
+  x'  cosa  -  y'  sma 


(94) 


6  -  «  -  (x'kcosak  -  y’ksinak)6  -  sinak6x, 

“c  yo 


-  cosa  6  , 
yc 


k*  =  vk  -  vk  +  x'ksinak  +  y:kcosak 


(95) 


y0  •  ^c 


[Pk  sinak  +  cosak]6q  +  [Qk  sinak  -  cosak]6p 
+  [cosak  (1  +  PkQk)  +  sinak  (Pk  -  Qk) ] 6 a 
=  -  [sinak  (1  +  PkQk)  -  cosak  (Pk  -  Qk) J 


(96) 


Other  variational  relationships  which  must  be  solved  simulta¬ 
neously  with  equations  (91)  through  (96)  (See  Appendix  B)  are: 


-102- 


1 


(a)  Variation  of  ligament  insertion  coordinates, 


i  -  6  +  (x\  sincxk  +  y  ^  cosak)6 

xj  o  J  J 

=  -  xk  +  xk  +  x^  cosctk  -  yj  sinak  (j=l,..,4) 


(97) 


V  k 

6  -  S  -  (xt  cosa  -  y^  sina  )6 

Yj  Y0  3  3  "a 

=  -  yk  +  yk  +  xj  sinak  +  yj  cosak  (j=l,..,4) 

(b)  Variation  of  articulating  surfaces,  yc  =  f^ (xc)  and 
_____ 
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[A2  -  2A3xk  ♦  3A4(xk)2  ♦  4A5Cxk)3]6Xc  -  Aj 
+  A2xk  +  A3(xk)2  +  A4(xk)3  +  A,(xk)4  -  yk 


lS^c' 


(99) 


V 

7  c 


(A£  +  2A^x^.k)6x,  =  Aj  +  A£x£k  +  A£x£k2  -  y£k 
c 


(100) 


(c)  Variation  of  first  derivatives,  P  and  Q  of  the  articulating 
surfaces,  y  =  f^(x)  and  y'  =  f2(x'),  respectively: 


6p  -  (2A3  +  6A4xk  +  12A5xk  )  <5 x 


+  4ArXk  -  Pk 
5  C 


A2  *  2AjX*  ♦  3A4x^2 


(101) 


60  '  2A3  6x»  =  A2  +  2a3  x’k  -  Qk 


(102) 
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(d)  Variation  of  the  components  of  unit  normals,  n^  and  n^: 


Cl  +  Pk  )3/2 


5p  = 


+  pk%l/2 


(1  +  ) 


(103) 


+  Pk  V/2 


6P  " 


k2  1/2  ’ 

(1  +  PK  )i/z 


g _  *  =  -gQ  .  J 

k2  3/2  Q  k2  1/2  2x' 

(1  +  Qk  )i/ZJ  (1  +  Qk  )i/Z 


(104) 


(105) 


. eg _  6  =  s  .  nk 

k2  3/2  Q  k2  1/2  2?' 

(1  +  Qk  )3/z  (1  +  Qk  )1/z 


(106) 


where  3  =  — rr 
dx  L 


d2f2  .  d2f! 
— ry  ,  and  y  =  — rj 

dx  z  dx  z 


Equations  (91)  through  (106)  form  a  set  of  22  simultaneous 
algebraic  equations  defined  by  equation  (52).  The  vector.  A, 
of  equation  (52)  has  the  following  elements:  SXq,  &yQ,  Sa>  <$N» 
5xc>  6Xj  (j  =  l,..,4),  Sy-j  ( j  =  1 ,  •  .  ,4)  ,  <Snlx,  6nly,  6n2xi  > 

<sn2yi  »  6yc,  5y^ >  6p  and  6q,  where 

<$x  =  x-coordinate  variation  of  the  center  of  mass  of  the 

0  tibia. 

<5  =  y-coordinate  variation  of  the  center  of  mass  of  the 

y  o  .  ... 

tibia. 

Sa  =  variation  of  the  orientation  angle,  a,  of  the  tibia 
with  respect  to  the  femur. 

=  variation  of  the  normal  force,  N. 

<5X  =  x-coordinate  variation  of  the  contact  point,  C,  in  the 

(x,y)  system. 

6X,  =  x-coordinate  variation  of  the  contact  point,  C,  in  the 

(x^y1)  system. 
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=  x-coordinate  variation  of  the  ligament  insertion 
points ,  ( j  =  1 , . . 4) . 

=  y-coordinate  variation  of  the  ligament  insertion 
points ,  ( j=l , . . 4) . 

=  x-component  variation  of  the  unit  normal  vector  to 
the  articulating  surface  of  the  femur. 

=  y-component  variation  of  the  unit  normal  vector  to 
the  articulating  surface  of  the  femur. 

=  x-component  variation  of  the  unit  normal  vector  to 
the  articulating  surface  of  the  tibia. 

=  y-component  variation  of  the  unit  normal  vector  to 
the  articulating  surface  of  the  tibia. 

=  y-coordinate  variation  of  the  contact  point,  C,  in 
the  fixed  (x,y)  system. 

=  y-coordinate  variation  of  the  contact  point,  C,  in 
the  moving  (x'jy1)  system. 

=  variation  of  the  first  derivative  of  the  femoral 
articulating  surface  equation. 

=  variation  of  the  first  derivative  of  the  tibial 
articulating  surface  equation. 

After  applying  an  external  force  and  moment  to  the  center 


of  mass  of  the  tibia,  the  above  system  of  22  equations  are 


solved  for  the  22  variables  by  the  computer  program  JNTMDL. 


SOME  ASPECTS  OF  THE  COMPUTER  PROGRAM,  JNTMDL 

Following  the  numerical  procedure  described  previously, 
the  system  of  22x22  equations  defined  by  equation  (52)  are 
solved  using  the  JNTMDL  program  given  in  Appendix  C. 
Prescribed  constant  parameters  used  in  the  program  are: 

(a)  the  coordinates  (xj,yj)  and  (xt,yj)  of, 
respectively,  the  insertions  and  origins  of  the 
ligaments . 

(b)  stiffness  values  of  the  ligaments,  kj . 

(c)  mass,  M,  and  mass  moment  of  inertia,  I  ,  of  the 
moving  tibia. 
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(d)  initial  linear  (xQ,  yQ ,  xQ ,  yQ)  and  angular  (a, 
a)  velocities  and  accelerations  of  the  center  of 
mass  of  the  moving  tibia. 

(e)  coefficients,  and  A^,  of  the  articulating 
surface  equations,  y  =  f^(x)  and  y'  =  f2(x'), 
respectively. 

(f)  applied  external  force,  F_ ,  and  moment,  M  . 

(g)  time  increment,  At. 

(h)  convergence  criterion,  omega. 

By  specifying  the  initial  contact  points,  xc  and  x'c ,  the 
JNTMDL  program  first  determines  the  starting  configuration  of 
the  moving  tibia  relative  to  the  fixed  femur.  This  is  done  by 
satisfying  the  geometric  compatibility  equations  (14)  and 
contact  conditions  (36)  .  After  applying  the  prescribed 
external  force  and  moment,  by  the  use  of  iteration  process, 
the  JNTMDL  program  calculates  the  22  delta  variations  at  each 
iteration.  The  iteration  process  at  a  fixed  time  station 
continues  until  the  delta  quantities  of  all  variables  become 
negligibly  small.  In  the  present  work,  a  solution  is  accepted 
and  the  iteration  process  is  terminated  when  the  delta 
quantities  become  less  than  or  equal  0.01%  of  the  previous 
values  of  the  corresponding  variables.  The  converged  solution 
of  each  variable  is  then  used  as  the  initial  value  for  the 
next  time  step  and  the  process  is  repeated  for  consecutive 
time  steps.  In. our  application  of  the  method,  only  5-6  itera¬ 
tions  were  necessary  for  convergence  most  of  the  time.  But 
there  were  also  cases  where  more  iterations  were  required  for 
convergence.  This  was  especially  true  at  instants  at  which 
there  was  a  sudden,  sharp  change  in  the  response  of  the  tibia. 
Such  behavior  manifested  itself  usually  when  the  tibia  started 
mqving  in  the  opposite  direction  due  to  the  fact  that  the 
pulling  force  of  a  ligament(s)  became  the  governing  force  of 
the  problem.  Shorter  time  steps  required  fewer  iterations, 
as  expected,  even  at  the  situation  described.  The  time 
increment  used  in  the  present  work  is  At  =  0.0001  second. 
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At  each  converged  solution,  the  following  quantities  are 
determined  by  JNTMDL  program: 

(a)  The  orientation  angle,  a  (flexion  angle),  of  the 
moving  tibia  relative  to  the  fixed  femur. 

(b)  The  position  (xQ,y0)  of  the  center  of  mass  of 
tibia. 

(c)  The  linear  (xQ , yQ)  and  angular  accelerations  of 
the  center  of  mass  of  the  moving  tibia. 

(d)  The  location  of  contact  points,  xc  and  x^.,  on  the 
fixed  femur  and  moving  tibia,  respectively. 

(e)  The  magnitude  of  the  contact  force,  N,  exerted  on 
the  moving  tibia. 

(f)  The  magnitude  of  ligament  forces,  F ^ . 

(g)  The  external  energy  supplied  to  and  internal  energy 
generated  by  the  moving  tibia,  for  comparison 
purposes . 

The  final  section  of  this  report  will  present  results 
from  this  program  for  several  loading  conditions  applied 
through  the  center  of  mass  of  the  tibia. 

NUMERICAL  RESULTS 

The  numerical  results  to  be  presented  are  only  for  an 
external  force  acting  on  the  tibia  without  the  presence  of  an 
external  moment.  It  is  assumed  that  the  force  is  always 
perpendicular  to  the  longitudinal  axis  of  the  tibia  (y'-axis) 
and  passes  through  its  center  of  mass.  Let  this  force  be 
denoted  by  P(t) .  A  parametric  study  of  the  effect  of  various 
combinations  of  moment  and  force  acting  simultaneously  on  the 
response  of  the  knee  joint  may  prove  to  be  rewarding. 

However,  for  the  present  work  we  will  only  consider  an  exter¬ 
nal  force  and  believe  that  this  will  be  sufficient  to 
illustrate  the  capabilities  of  the  model.  The  effect  of  the 
shape  of  the  forcing  function  on  the  knee  joint  response  will 
be  studied  by  considering  the  following  two  functions  for  P(t): 
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(107) 


P(t)  -  A[H(t)  -  H(t-t0)] 

which  is  a  rectangular  pulse  of  duration  tQ  ,  and  amplitude  A; 
and 


-4.73(t/tflr  /t\ 

P(t)  =  Ae  u  sin  (108) 

which  is  an  exponentially  decaying  sinusoidal  pulse  of 
duration  t^,  and  amplitude  A.  A  dynamic  loading  in  the  form 
of  equation  (107)  is  extremely  difficult  to  simulate 
experimentally;  however,  the  study  of  these  two  functions  will 
hopefully  be  helpful  in  understanding  the  effect  of  rise  time 
of  the  dynamic  load  on  the  joint  response.  Equation  (108)  is 
a  more  realistic  forcing  function  and  it  has  been  previously 
used  as  a  typical  representation  of  the  dynamic  load  in  head 
impact  analysis  (Engin  and  Akkas ,  [1978]). 

The  following  are  obtained  as  a  function  of  time  from 
the  computer  program,  JNTMDL;  the  coordinates  xQ,yo  of  the 
center  of  mass  of  tibia;  the  flexion  angle  a;  the  coordinates 
xc  and  x^  of  the  contact  point  in  (x,y)  and  (x',y')  coordinate 
systems  respectively;  the  magnitude,  N,  of  the  contact  force; 
the  elongations  of  the  ligaments  and  the  ligament  forces,  Fj ; 
and  the  internal  and  external  energies  of  the  system. 

The  initial  values  for  xQ,  yQ  and  a  are  obtained  by 
specifying  the  location  of  the  starting  contact  point.  Here, 
the  following  values  are  used  for  the  coordinates  of  the 
contact  point  at  t  =  0: 

xc  =  -4.2  cm  ,  x^  =  2.5  cm 
which  yields 

xQ  =  -20.16  cm  ,  yQ  =  17.49  cm  ,  and  a  =  234.79° 

This  angle  of  rotation  a,  corresponds  to  a  flexion  angle  of 
54.79°  for  which  the  ligaments  of  the  knee  joint  are  in  a 


relatively  relaxed  condition.  It  was  reported  as  early  as 
1907  by  Pringle  that  the  position  of  maximum  relaxation  of  the 
knee  joint  ligaments  was  approximately  halfway  between  full 
flexion  and  full  extension  of  the  knee  joint. 

The  effect  of  pulse  duration  on  the  response  of  the  knee 
joint  motion  is  studied  by  taking  t  =  0.05S,  0.10S,  and  0.15S 
for  both  rectangular  and  exponentially  decaying  sinusoidal 
pulses.  The  effect  of  pulse  amplitude,  is  also  examined  by 
taking  A  =  20N,  60N,  100N,  140N  and  180N  for  both  types  of 
pulses . 

Ligament  forces  as  functions  of  flexion  angle  of  the  knee 
joint  for  the  two  previously  described  forcing  functions  are 
presented  in  Figures  35  through  40.  Results  indicate  that 
when  the  knee  joint  is  extended,  by  a  dynamic  application  of 
a  pulse  on  the  tibia,  lateral  collateral,  medial  collateral 
and  anterior  cruciate  ligaments  are  elongated  while  the 
posterior  cruciate  ligament  is  shortened.  The  load  carried  by 
the  anterior  cruciate  ligament  is  substantially  higher  than 
those  of  the  lateral  collateral  and  medial  collateral 
ligaments . 

The  variation  of  the  lengths  of  ligaments  and  the  forces 
carried  by  them  during  normal  knee  motion  has  been  the  subject 
of  various  studies  reported  in  the  literature  and  in  these 
studies  several  different  opinions  and  conclusions  have  been 
expressed  in  regard  to  the  biomechanical  role  and  function  of 
various  ligaments  of  the  knee  joint.  The  function  of  the 
anterior  cruciate  as  depicted  in  the  dynamic  knee-model 
developed  in  this  research  program  is  to  resist  anterior  dis¬ 
placement  of  the  tibia.  This  function  is  in  general  agreement 
with  the  experimental  and  clinical  studies  of  Kennedy  and 
Fowler  [1971];  Girgis,  Marshall  and  Monajem  [1975];  Van  Dijk, 
Huiskes  and  Selvik  [1979];  and  quasi-static  model  analyses  of 
Crowninshield,  Pope  and  Johnson  [1976];  and  Wismans  [1980]. 
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Figure  35.  Ligament  forces  as  functions  of  flexion  angle, 
for  an  externally  applied  rectangular  pulse  of  0.05  second 
duration  and  amplitude  of  (a)  60  N,  (b)  100  N,  (c)  140  N 
and  (d)  180  N. 
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Figure  36.  Ligament  forces  as  functions  of  flexion  angle, 
for  an  externally  applied  rectangular  pulse  of  0.10  second 
duration  and  amplitude  of  (a)  60  N,  (b)  100  N,  (c)  140  N 
and  (d)  180  N. 
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Figure  37.  Ligament  forces  as  functions  of 
for  an  externally  applied  rectangular  pulse 
duration  and  amplitude  of  (a)  60  N,  (b)  100 
and  (d)  180  N. 


flexion  angle, 
of  0.15  second 
N,  (c)  140  N 
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Figure  38.  Ligament  forces  as  functions  of  flexion  angle, 
for  an  externally  applied,  exponentially  decaying  sinusoidal 
pulse  of  0.05  second  duration  and  amplitude  of  (a)  60  N, 

(b)  100  N,  (c)  140  N  and  (d)  180  N. 
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Figure  39.  Ligament  forces  as  functions  of  flexion  angle, 
for  an  externally  applied,  exponentially  decaying  sinusoidal- 
pulse  of  0.10  second  duration  and  amplitude  of  (a)  60  N, 

(b)  100  N,  (c)  140  N  and  (d)  180  N. 
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Figure  40.  Ligament  forces  as  functions  of  flexion  ai\gle, 
for  an  externally  applied,  exponentially  decaying  sinusoidal 
pulse  of  0.15  second  duration  and  amplitude  of  (a)  60  N, 

(b)  100  N,  (c)  140  N  and  (d)  180  N. 


-115- 


The  role  of  the  posterior  cruciate  as  predicted  by  the 
model  is  to  resist  posterior  displacement  of  the  tibia.  This 
function  of  the  posterior  cruciate  ligament  is  in  agreement 
with  the  experimental  studies  of  Kennedy  and  Grainger  [1967] ; 
Edwards,  Lafferty  and  Lange  [1970];  Girgis,  Marshall  and 
Monajem  [1975]  ;  Crowninshield,  Pope  and  Johnson  [1976]  and 
Wismans  [1980] . 

The  present  dynamic  model  also  predicts  that  the  medial 
and  lateral  collateral  ligaments  offer  very  little  resistance 
in  the  flexion-extension  motion  of  the  knee  joint.  The  major 
role  of  these  ligaments  is  to  offer  varus-valgus  and  partial 
internal-external  rotational  stability.  The  model  shows  that 
as  the  knee  joint  is  extended  under  influence  of  a  dynamic 
load  the  lateral  collateral  and  medial  collateral  ligaments 
elongate  at  different  magnitudes.  This  prediction  is  in 
general  agreement  with  the  experimental  results  of  Smillie 
[1970] ;  Edwards,  Lafferty  and  Lange  [1970]  ;  and  Wang,  Walker 
and  Wolf  [1973] . 

The  model  shows  good  agreement  with  the  quasi-static 
experimental  investigations  reported  in  the  literature.  It 
is  important  to  note  that  the  dynamic  model  presented  in  this 
report  is  an  idealized  representation  of  a  very  complex 
anatomical  structure;  thus,  static  experimental  studies  may 
not  support  some  of  the  predictions  of  the  model.  Additional 
disagreements  may  also  be  due  to  approximate  locations  of  the 
attachment  sites  of  the  ligaments  in  particular,  and  two- 
dimensional  nature  of  the  model  in  general. 

In  Figures  41  and  42,  a  few  representative  plots  of 
forces  in  the  anterior  cruciate  and  lateral  collateral  liga¬ 
ments  are  plotted  as  a  function  of  time  for  two  different 
forcing  functions  with  varying  pulse  durations.  Although  not 
presented  here,  similar  curves  may  be  obtained  for  other  pulse 
durations  and  pulse  magnitudes.  Generally,  the  shorter  the 
pulse  duration,  for  a  fixed  amplitude,  the  sooner  the  tibia 
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(a)  (b) 

Figure  41.  Anterior  cruciate  ligament  force  as  a  function 
of  time  for  externally  applied  (a)  rectangular  and  (b)  ex¬ 
ponentially  decaying  sinusoidal  pulses,  of  60  N  amplitude 
and  durations  of  0.05,  0.10  and  0.15  second. 


(a) 


(b) 


Figure  42.  Lateral  collateral  ligament  force  as  a  function 
of  time  for  externally  applied  (b)  rectangular  and  (a)  ex¬ 
ponentially  decaying  sinusoidal  pulses,  of  60  N  amplitude 
and  durations  of  0.05,  0.10  and  0.15  second. 
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reaches  its  turning  point  (i.e.,  direction  of  motion  reverses) 
and  for  a  given  pulse  duration,  the  smaller  the  amplitude,  the 
sooner  the  turning  point  is  reached.  In  Figures  44-47  the 
values  in  parentheses  indicate  the  flexion  angles  at  the 
corresponding  times.  Note  that,  for  illustrative  purposes  up 
to  6°  of  hyperextension  was  allowed.  Generally,  one  expects 
only  1  to  3°  of  hyperextension  to  be  anatomically  tolerable 
beyond  which  joint  failure  becomes  unavoidable. 

In  Figures  43  and  44,  contact  forces  as  a  function  of 
time  are  plotted.  These  forces  are  in  response  to  the 
different  forcing  functions  with  varying  amplitudes  and  pulse 
durations.  Note  that  the  magnitudes  of  the  ligament  and  the 
corresponding  contact  forces  in  response  to  a  particular 
forcing  function  are  comparable. 

Femoral  and  tibial  contact  point  locations  as  a  function 
of  flexion  angle  are  plotted  in  Figures  45  through  50.  In 
these  figures  the  values  in  the  parentheses  indicate  the  total 
elapsed  time  of  the  motion  since  its  initiation.  As  the 
flexion  angle  decreases,  it  can  be  seen  that,  the  curves 
representing  femoral  contact  points  have  a  steadily  increasing 
positive  slope,  while  the  curves  of  tibial  contact  points 
change  slope  from  positive  to  negative  or  vice-versa  at 
various  flexion  positions.  This  phenomenon  may  be  explained 
by  the  combined  rolling  and  sliding  motion  of  the  tibia  on 
the  femur;  that  is,  positive  and  zero  slopes  of  the  curves 
representing  the  tibial  contact  points  can  be  interpreted  as 
corresponding  to  the  sliding  motion,  while  negative  slopes 
indicating  the  rolling  motion  of  tibia  on  femur.  Generally, 
the  curves  of  tibial  contact  points  have  predominately 
negative  slopes  toward  the  end  of  the  extension  motion 
indicating  that  in  this  part  of  the  motion  rolling  is  the 
'essential  component.  These  results  are  in  general  agreement 
with  the  work  of  Walker  and  Hajek  [1972];  Kettelkamp  and 
Jacobs  [1972];  and  Wismans,  et  al.  [1980]. 
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Figure  43.  Contact  forces  as  a  function  of  time  for  an 
externally  applied  rectangular  pulse:  (a)  60  N  amplitude 
and  durations  of  0.05,  0.10  and  0.15  second;  20  N,  60  N, 

100  N,  140  N  and  180  N  amplitudes  and  durations  of  (b)  0.05 
second;  (c)  0.10  second  and  (d)  0.15  second. 
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Figure  45.  Femoral  and  tibial  contact  points  as  a  function 
of  flexion  angle  for  an  externally  applied  rectangular  pulse 
of  0.05  second  duration  and  amplitudes  of  (a)  60  N,  (b)  100  N 
(c)  140  N  and  (d)  180  N. 
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Figure  46.  Femoral  and  tibial  contact  points  as  a  function 
of  flexion  angle  for  an  externally  applied  rectangular  pulse 
of  0.10  second  duration  and  amplitudes  of  (a)  60  N,  (b)  100  N 
(c)  140  N  and  (d)  180  N. 
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Figure  47.  Femoral  and  tibial  contact  points  as  a  function 
of  flexion  angle  for  an  externally  applied  rectangular  pulse 
of  0.15  second  duration  and  amplitudes  of  (a)  60  N,  (b)'  100  N 
(c)  140  N  and  (d)  180  N. 
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Figure  48.  Femoral  and  tibial  contact  points  as  a  function 
of  flexion  angle  for  an  externally  applied  exponentially 
decaying  sinusoidal  pulse  of  0.05  second  duration  and  ampli¬ 
tudes  of  (a)  60  N,  (b)  100  N,  (c)  140  N  and  (d)  180  N. 
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(a) 


(b) 


(c) 


(d) 


Figure  50.  Femoral  and  tibial  contact  points  as  a  function 
of  flexion  angle  for  an  externally  applied,  exponentially 
decaying  sinusoidal  pulse  of  0.15  second  duration  and  ampli¬ 
tudes  of  (a)  60  N,  (b)  100  N,  (c)  140  N  and  (d)  180  N. 
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In  Figures  51  and  52,  total  internal  energy  of  the  moving 
tibia  as  a  function  of  flexion  angle  are  plotted.  This  total 
internal  energy  is  defined  as  the  summation  of  potential 
energies  of  the  nonlinear  elastic  springs,  simulating  the 
ligament  forces,  and  the  kinetic  energy  of  the  tibia.  As  it 
can  be  seen  from  these  curves,  for  a  given  pulse  shape  of 
fixed  duration,  as  the  pulse  magnitude  increases,  the  knee 
joint  further  extended  and  the  entire  motion  has  a  shorter 
response  time  period.  External  energy  of  the  system  is  also 
calculated  based  on  the  applied  external  forces.  Continuous 
plots  of  total  internal  and  external  energy  of  the  system  for 
two  different  pulse  shapes  of  specified  magnitude  and  duration 
are  presented  in  Figures  53  and  54.  These  figures  show  that 
the  total  internal  and  external  energies  of  the  system  remain 
the  same  except  for  the  latter  part  of  the  motion,  which  is 
due  to  the  accumulation  of  numerical  round-off  errors. 

Finally,  for  illustrative  purposes,  with  the  aid  of  Versatec 
plotter,  continuous  plots  of  xQ,y0  (coordinates  of  the  center 
of  mass  of  the  tibia) ;  xc  and  (femoral  and  tibial  contact 
points,  respectively)  and  flexion  angle,  a,  as  a  function  of 
time  are  plotted  in  figures  55  and  56. 

SUMMARY  AND  CONCLUDING  REMARKS 

The  research  work  discussed  and  presented  in  this  report 
can  be  summarized  in  the  following  paragraphs: 

1.  Detailed  anatomical  description  and  articulation  of 
the  elbow,  shoulder,  hip,  knee  and  ankle  joints  with  suffi¬ 
cient  illustrative  figures  for  each  joint  were  presented  and 
major  anatomical  parts  have  been  identified  using  the 
generally  applied  medical  terminology. 

2.  Mathematical  descriptions  of  the  articulating 
surfaces  of  elbow,  hip,  knee  and  ankle  joints  have  been 
determined  by  means  of  a  sonic  digitizing  technique.  The 
attachment  sites  of  the  major  ligaments  of  each  joint  were 
also  determined  ‘and  tabulated. 
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Figure  51.  Total  internal  energy  of  the  moving  tibia  as  a 
function  of  flexion  angle  for  an  externally  applied  (a)  rec 
tangular  pulse  and  (b)  exponentially  decaying  sinusoidal 
pulse  of  0.10  second  duration  and  amplitudes  of  20  N,  60  N, 
140  N  and  180  N. 


(a)  (b) 

Figure  52.  Total  internal  energy  of  the  moving  tibia  as  a 
function  of  flexion  angle  for  an  externally  applied  (a)  rec 
tangular  pulse  and  (b)  exponentially  decaying  sinusoidal 
pulse  of  0.15  second  duration  and  amplitudes  of  20  N,  60  N, 
100  N,  140  N  and  180  N. 
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Figure  53.  Continuous  total  internal  and  external  energies 
of  the  joint  system  as  a  function  of  time  for  an  externally 
applied  rectangular  pulse  of  20  N  amplitude  and  0.15 
second  duration 


Figure  54.  Continuous  total  internal  and  external  energies 
of  the  joint  system  as  a  function  of  time  for  an  externally 
applied,  exponentially  decaying  sinusoidal  pulse  of  100  N 
amplitude  and  of  0.15  second  duration. 


Figure  55.  Continuous  plots  of  tibial  center  of  mass  co¬ 
ordinates  (XO,YO) ,  femoral  (XC)  and  tibial  (XPC)  contact 
points  and  the  flexion  angle  (ALFA)  as  functions  of  time 
for  an  externally  applied  rectangular  pulse  of  20  N  am¬ 
plitude  and  of  0.15  second  duration 
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Figure  56.  Continuous  plots  of  tibial  center  of  mass  co¬ 
ordinates  (XO,YO),  femoral  (XC)  and  tibial  (XPC)  contact 
points  and  the  flexion  angle  (ALFA)  as  functions  of  time 
for  an  externally  applied,  exponentially  decaying  sinusoidal 
pulse  of  100  N  amplitude  and  of  0.15  second  duration. 


3.  A  review  of  the  available  literature  on  the 
biomechanical  behavior  of  soft  tissues  in  general  and  liga¬ 
ments  in  particular  was  presented.  An  appropriate  constitu¬ 
tive  equation  for  the  elastic  behavior  of  the  ligaments  was 
established. 

4.  Two-  and  three-dimensional  mathematical  dynamic 
models  of  a  general  two-body-segmented  articulating  joint 
have  been  formulated  in  order  to  describe  the  relative  motions 
between  the  segments  and  the  various  forces  produced  at  the 
joint.  The  governing  equations  for  these  models  are  set  of 
highly  nonlinear  equations  and  numerical  solutions  were 
discussed  in  some  detail. 

5.  The  two-dimensional  model  was  applied  to  the  knee 
joint  and  the  numerical  results  from  this  model  were  presented 
to  illustrate  the  effects  of  duration  and  shape  of  the 
dynamically  applied  loads  on  the  response  of  the  joint. 

Special  attention  has  been  given  to  the  ligament  and  contact 
forces,  the  location  of  contact  points,  anterior-posterior 
displacements  and  the  comparison  between  the  internal  and  the 
external  energy  of  the  system.  The  results  were  compared 
with  the  available  experimental  data  from  the  literature  to 
establish  the  validity  of  the  model. 

It  is  appropriate  to  make  several  remarks  on  the 
numerical  techniques  tried  in  the  course  of  obtaining  the 
solution  for  the  governing  equations  of  the  two-dimensional 
dynamic  model  (equations  (45),  (46)  and  (47),  coupled  with 
constraint  equations  (14)  and  (36)).  In  the  first  method,  the 
second-order  differential  equations  were  transformed  to  a  set 
of  nonlinear  algebraic  equations  by  substituting  for  the 
differential  elements,  their  equivalent  backward  difference 
approximations.  In  this  case,  it  was  impossible  to  obtain  a 
converging  solution  due  to  the  highly  nonlinear  nature  of 
these  equations. 
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In  the  second  method,  the  flexion  jangle,  a,  of  the  moving 
tibia  and  contact  point  coordinates,  xc  and  x'c  were  obtained 
from  the  simultaneous  solution  of  the  three  nonlinear  con¬ 
straint  equations.  Knowing  a,  a  was  obtained  via  backward 
difference  method  and  then  the  normal  force,  N,  was  determined 
from  the  third  equation  of  motion.  The  other  two  second- 
order  nonlinear  differential  equations  which  were  in  terms  of 
x„  and  v  were  written  as  a  set  of  four  first  order 
differential  equations  by  direct  substitutions.  Runge-Kutta 
method  was  applied  to  these  equations  and  solutions  for  xQ, 
yQ  and  their  time  derivatives  were  obtained.  Using  the  new 
values  of  xQ  and  y  ,  a  new  value  of  a  and  contact  point 
coordinates  were  obtained  and  the  entire  procedure  was  repeated 
for  the  next  time  step.  Although  mathematically  all  the 
geometric  constraints  and  governing  equations  of  motion  were 
satisfied,  the  results  obtained  using  this  second  method  were 
not  in  agr  ment  with  the  actual  physical  geometry  and  the 
anatomy  of  the  joint.  This  was  concluded  to  be  due  to  the 
solution  technique  which  was  not  solving  all  the  equations 
simultaneously,  and  during  the  solution  process  it  was  forcing 
some  of  the  variables  to  accept  values  which  were  mathemati¬ 
cally  correct  but  physically  unacceptable.  Finally,  the 
Newton-Raphson  iteration  process  along  with  Newmark  method  of 
differential  approximation  was  chosen  as  the  method  of 
solution  which  yielded  accurate  and  stable  solutions  for  the 
model.  The  entire  numerical  procedure  has  been  explained  in 
detail  in  previous  sections. 

The  extensions  of  the  research  work  presented  in  this 
report  can  be  in  several  areas.  One  can  investigate  the 
influence  of  the  variations  of  initial  strain  of  the  ligaments 
and  their  attachment  sites  (i.e.  insertions  and  origins)  on 
the  response  of  the  model.  A  parametric  study  addressing  to 
these  points  may  reveal  the  sensitivity  of  the  model  to  the 
variations  of  the  coordinates  of  the  insertions  and  origins  of 
the  ligaments. 
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In  a  similar  way  the  mathematical  models  of  the  other 
major  articulating  joints  can  be  developed.  However,  a 
special  attention  should  be  given  in  modeling  of  the  ligaments. 
Some  ligaments,  particularly  the  thick  band  or  large  cord-like 
ligaments,  have  complex  behavior,  with  various  portions 
behaving  differently  under  given  conditions  or  configurations. 
These  are  described  in  the  literature  as,  in  the  case  of  the 
broad  ligaments  of  the  hip  joint,  having  anterior  and 
posterior  fibers,  or  medial  and  lateral  components.  The 
mathematical  model  should  include  additional  elastic  elements 
to  reflect  the  contributions  of  various  fibers  of  the  ligaments. 
Unfortunately  proper  experimental  data  to  determine  the 
constitutive  behavior  of  these  thick  band  or  broad  ligaments 
do  not  exist. 

Finally,  the  numerical  procedure  utilized  in  the  solution 
of  the  two-dimensional  dynamic  joint  model  equations  should 
be  applied  to  the  three-dimensional  dynamic  model  formulation 
presented  in  this  report.  Considering  importance  and  rele¬ 
vance  of  a  three-dimensional  dynamic  joint  model  to  the  Air 
Force  related  applications,  a  very  serious  effort  should  be 
expended  in  the  direction  of  obtaining  numerical  solutions 
for  this  complex  joint  model. 
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APPENDIX  A 


COMPUTER  SUBROUTINE  CHEPLS: 

A  LEAST  SQUARES  CURVE-FITTING  ROUTINE 

The  calculation  most  commonly  performed  on  experimen¬ 
tal  data  is  to  fit  the  data  with  a  polynomial  of  the  form: 
F(X)  =  A  +  BX  +  CX2  +  DX3  +  ...  MXN,  which  is  best  in  the 
least  squares  sense.  CHEPLS  routinely  performs  this  cal¬ 
culation.  Data  input  is  simple  enough  so  that  no  prior  com¬ 
puter  experience  is  required  for  its  use. 

The  routine  determines  by  means  of  statistical  tests 
whether  the  set  of  data  is  linear  or  non-linear.  If  the 
data  are  non-linear  at  the  95  percent  confidence  level,  the 
routine  finds  the  lowest  degree  polynomial  which  adequately 
represents  the  data.  The  calculation  procedure  is  to  com¬ 
pare  the  standard  deviation  of  the  model  of  degree  n,  with 

a.  1- 

the  standard  deviation  of  the  (n+1)  degree  model.  If 
there  is  no  significant  difference  at  the  95  percent  con¬ 
fidence  level,  then  it  can  be  said  that  in  19  out  of  20 
such  sets  of  data,  the  polynomial  of  degree  n  is  the  low¬ 
est  degree  polynomial  which  adequately  represents  the  data. 

CHEPLS  prints  out  the  least  squares  coefficients  for 
all  polynomials,  from  the  first  degree  through  the  highest 
degree  calculated.  For  each  polynomial,  the  standard  de¬ 
viation  and  the  confidence  bands  about  each  coefficient  at 
the  95  percent  level  are  printed.  Values  of  Y  are  calcu¬ 
lated  from  the  model  for  every  value  of  X,  and  are  compared 
with  the  observed  (ie,  experimental  or  input)  values  of  Y. 
The  least  squares  matrix  is  printed  for  reference  in  the 
event  additional  statistical  calculations  are  needed.  The 
Y  values  are  coded  by  subtracting  a  constant  to  reduce 
round-off  error.  As  many  as  1000  data  points  may  be  sub¬ 
mitted  in  each  data  set.  CHEPLS  will  produce  polynomials 
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up  to  degree  ten.  Data  sets  may  be  stacked,  one  behind  the 
other  for  additional  data  sets. 

RESTRICTIONS 

(a)  The  number  of  data  points  must  be  less  than  or  equal  to 

1000. 

(b)  The  maximum  degree  polynomial  is  10.  Even  in  double¬ 
precision,  round-off  error  may  accumulate  sufficiently 
to  invalidate  results  for  even  as  low  as  a  5t^1  degree 
polynomial . 

(c)  All  numeric  data  must  be  right- justified  in  the  field. 
INPUT  DATA 

Each  data  set  consists  of  a  data  header,  a  label  and 
one  or  more  data  cards  and  a  termination  card.  For  stacked 
data  sets,  only  the  final  data  set  may  contain  a  termination 
card.  All  fixed  point  variables  or  integers  are  entered 
without  a  decimal.  Floating  point  variables  require  a  deci¬ 
mal  in  the  field. 

(a)  Data  Header 
Columns  4-5 

Columns  9-10 


(b)  Label 

Columns  1-80 


Enter  the  number  of  data  points 
in  this  set. 

Blank  for  normal  use.  For  ob¬ 
taining  all  polynomials  up  to  the 
n*^'  degree  (maximum  of  10) ,  enter 
n.  All  least  squares  polynomials 
will  then  be  calculated  for  degrees 
1  through  n. 

Enter  any  alphanumeric  identifica¬ 
tion  to  be  included  on  the  print¬ 
out  for  descriptive  purposes. 


4 
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(c)  Data  Cards 

Columns  1-80  Enter  two  data  pairs  per  card,  as 

X,  Y,  X,  Y.  Format  may  be  either 
4E18.8  or  4F18.8.  Note  that  for 
an  odd  number  of  data  sets,  columns 
37-80  of  the  final  data  card  will 
be  blank. 

(d)  Termination  Card 

Columns  23-25  Enter  END  to  signify  the  end  of  da¬ 

ta  entry:  no  further  data  sets. 

The  job  control  language  for  running  this  subroutine  is  as 
follows : 


//  TIHE=(0>25) »REGI0N=192K 

/BJOBPARM  L I NES= 1 000 » CARDS9 1 00 » DISK I0« 1300 

//« 

//*  EXECUTE  THE  CHEH  ENG  LIBRARY  PROGRAM  CHEPLS 

//* 

Y/STEPB  EXEC  PGM=IEHL»PARM='XREF,LIST»MAP'rTINE*(0»30) 

//*  STEPB  IS  LINK  EDIT  STEP  (LOAD  NODULE  ASSEMBLED) 

//SYSLIB  DD  DSNAME=SYS1.F0RTLIB.DISP*SHR 
//  DD  DSNAME*FEAA80.CHEMENGR»DISP*SHR 

//  DD  DSNAME*SYS2< F0RTSSP » DISP=SHR 

//0LDLIB  DD  DSNAME*FEA680 • CHEMENGR » DISP=SHR 

//SYSLMOD  DD  DSNAME*I100(MAIN)»UNIT»SYSDAiSPACE*(CYL>(l»l»l))»  * 

//  D I SP» ( NEW » PASS ) » DCB* ( RECFM»U » BLKS I ZE-3072 ) 

//SYSPRINT  DD  SYS0UT-A 

//SYSUT1  DD  UNIT«SYSDA.SPACE»(CYL>(2»1)) 

/ /SYSLIN  DD  t 

INCLUDE  OLDLID(CHEPLS) 

/» 

//STEPC  EXEC  PGM-t. STEPB. SYSLM0D»RE6I0N-126K»TINE-(lf00) 

//«  STEPC  IS  EXECUTION  STEP 
//SYSUDUMP  DD  SYSOUT-A 
//FT0SF001  DD  DDNAME'SYSIN 
//ITMF001  DD  SYSOUT-A 
//FT07F001  DD  SYS0UT>B 
VSYSIN  DD  I 

ft*  REMOVE  THIS  CARD  AND  REPLACE  UITK  DATA  DECK  (BLANK  CARD  IF  NO  DATE) 
/t 

// 
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b)  From  equation  (46)  : 


7  F.  +  YN  n.  +  (F  ) 

jy  ly  e  y 


=  My, 


in  (B3)  let: 


N  =  N  +  6 


N 


Then 


"ly  =  nly  +  6nly 


y  =  y  +fi 

7  o  7  o  y 


J=1  Fj y  ♦  Y[(Nk  ♦  eN)(nkly  +  6nly)]  ♦  (Fe)y 
l  Fjy  +  y[Nknky  +  Nk6nlv  +  nkv  «N]  +  (FJ 


j-1 


’nly  “ly  UNJ 


Let 


°  =l 


y„  kyo  *  «y  )  -  C"!  -  it  V 

'  o 


At 


Then 


(mjy  )«N  *  (YNk)«nly  -  (-^V  -  -  I  F 


At  70  j  =1 


..k  k  j  4  r  k  t  -  A 1 1  4  •  t  -  £ 

yN  niy  * M  |^r  (y0  -  yo  i  -  if  y0 


c)  From  equation  (47) : 

i^j  -  V>Fjy  -  yo>Fjx  *  y"t(x( 


■(yc  '  V"!*1  *  Me  =  V 


CBS) 


Then 


!t'{: 


IC<*k  +  O 

At'2  a 


t-At  4  •  t- At  "t-At-i  l 

a  rt  °  '  a  i  J 


"  S./jy  *  TNknly1Sx0 

,k  k 


+  U/jx  +  yN  nlx^yo  +  ^xc  -  xXy  ‘  "  Mx^N 


,k,..k 


-  lTN“(y*  -  #]4  *  bNk(xk  -  x*)]Jnly  *  [YN*nJ)rJsx 
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CONSTRAINT  EQUATIONS 
From  equation  (8) 


x  =  x  +  x’  Cosct 
c  o  c 


-  y'c  Sina 


(B7) 


in  ( B 7 )  let 


x  =  x  +  6 

C  C  X 


X  =  X  +  6 
o  O  X 


xc  =  xc  +  V 

c 

yc  =  yc  +  V 
7  c 


a  =  a  +  6 


Then 


(x*  +  «  ) 

c  V 


(xo  +  6x  >  +  (Xrk  +  «v')Cos(ak  +  6  ) 
v->  a  ^  a  a 


C/^k  +  6yOSin(ak  +  «o) 


assuming 


Cos(ak  +  6  )  =  Cosak  -  6  Sinak 
a  a 

Sin(ak  +  6  )  =  Sinctk  -  6  Cosak 
a  a 


then 
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*-Xc  +  Sx  ^  =  (xo  +  6X  ^  +  Cxck+  6x')(CoSak  -  «  Sinak) 
c  o  c  a 

-  (y'k+  <5  .)(Sinak  +  6  Cosak) 

c  yc  a 

or 

(xo  +  5x  )  +  x^kCosak  -  x'kSinak6a  +  CoSak6x, 

0  xc 

-  y'ksinak  -  y'kCosak6  -  Sinak6  , 

c  c  a  yc 

6xc  '  6xq  +  (x^k  Sinak  +  y£k  Cosak)6a  -  Cosak6x, 

c 

+  Sinak5y,  =  xk  -  xk  +  x^k  Cosak  -  y£k  Sinak  (B8) 

Similarly  from  equations  (8): 

yc  =  yo  +  xc  Sina  +  yc  Cosa  (B9) 

the  variation  equation  is: 

6yc  "  6yo  ’  (XCkC0Sak  '  ycksin“k)6a  •  Sin“\. 

-  Cosak5  ,  =  yk  -  yk  +  x,kSinak  +  y'kCosak 
y^/o/c  c  7c 

From  equation  (36) 

•  /  d£l  df 7  \  /dfi  df7\ 

Slna +  33 r~  “cfirrJ  "  Cosa  yxir  *  J =  0 

From  eqs.  (85),  (86)  and  (87),  equation  (Bll)  can  be  written  as 

Sina [1  +  (A2  +  2A3xc  +  3A4x^  +  4A5xc3)(A£  +  A'x^)] 

-  Cosa[(A2  +  2  AjXc  +  3A4xc2  +  4A5xc3)-(A^  +  2A£  xp  ]  =  ’0  (B12) 


(BIO) 


(Bll) 
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in  (B12)  let: 


P  =  A-  +  2A_x  +  3A.x  2  +  4Acx  3 
2  3  c  4  c  5  C 

Q  =  A'  +  2Aix ' 

x  2  3  c 

therefore,  the  normal  constraint  equation  (B12)  becomes: 

Sina[l  +  PQ]  -  Cosa[P-Q]  =  0 
in  (B13)  let: 

P  =  I*  +  6p 

Q  -  Qk  +  «Q 

k  .  . 

a  =  a  +6 

a 

Then 

Sin(ak  +  6a)[l  +  (Pk  +  6p)(Qk  +  6Q)]-Cos(ak  ♦  6a)|Pk  + 
6p  -  Qk  -  «Q]  =  0 

or 

(Sinak  +  6aCosak)[l  +  P kQk  +  Pk6Q  +  Qk5pl  -  (Cosak  - 
6a  Sinak) [Pk  -  Qk  +  «p  -  «q]  “  0 

or 

[Sinak  +P  kQkSinak  +  Pk  Sinak6q  +  Qk  Sinok  6p  +  Cosak 
+  RkQk  Cosok6  ]  +  [ -CosakDk  +  k  SinakS  +  Qk  Cosak 
-  Qk  Sinak6Q  -  Cosak6p  +  Cosa^Sq]  =  0 


(B13) 


6 

a 


or 


[P"  Sina*  +  Cosa *]6q  +  [Qk  Sinak  -  Cosa*]^ 
+  [Cosak  +  PkQk  Cosak  +  Pk  Sinak  -  Qk  Sinctk]6 

a 

+  [Sinak  +  PkQk  Sinak  -  Pk  Cosak  +  Qk  Cosak]  = 


finally: 


[pk  Sinak  +  Cosak]  Sq  +  [Qk  Sinak  -  Cosak]6p 
+  [Cosak  (1  -  P  kQk)  +  Sinak  (Pk  -  Qk)]o 

a 

=  -  [Sinak  (1  +  PkQk)  -  Cosak  (Pk  -  Qk)  ] 


(B14) 


P  =  A-  +  2A,x„  +  3A.X  2  +  4Arx  3 


in  CB15)  let: 


'P  -  P  +  6, 


Xc  =  xc  +  6x 


Then 


CB15) 


P  +  «p  "  a2  +  2A3(xk  +  6x  )  ♦  3A4(xk  +  :x  )2  +  4A5(xk  +  «x  ) 


6p  -  (2A-  +  6A.xX  +  12ApXk  )« 


5  A.  _  j 

C  '  X. 


A2  *  2Ajxk  *  3A4x^' 


k"3  k 
+  4ACXX  -  PK 

5  C 


(B16) 


Similarly  for: 


Q  «  A’  +  2Aix 1 

L  J  C 
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3(1  +  Q  )  -  3Q  <5, 


(1  +  Qk  )3/2 


and  finally 


2yt  Ld  +  Qk  )3/2J  (1  +  Qk  )1/2 


LIGAMENT  ATTACHEMENT  COORDINATES 

From  equation  (8),  the  ligaments  insertion  coordinates  are: 


Xj  =  xQ  +  xj  Cosa  -  yj  Sino 


y.  =  y  +  x!  Sina  +  y!  Cosa 
3  o  j 


in  equation  (B26)  let 


then 


*  s 

Xk  +  6 

o  x_ 


k  *  * 
a  +5 


(xj  +  6X.}  =  (xo  +  6x0)  +  xj  Cosc«k  +  -  yj 

=  (xk  +5  )  +  x! (Cosak  -  6  sinak) 

v  o  x  '  J  v  a  J 


+  y ! (Sinak  +  6  Cosak) 


4 


6  -6  +  (x!  Sinak  +  y!  Cosak)6 

x .  x  v  l  '  ]  t 

j  o  J  J 


xk  +  xk  +  xj  Cosak  -  yj  Sinak  (j=l,..,4) 


(B28) 


similarly,  in  equation  (B27)  let: 


yk  +  <5 
J  y 

yk  +  6 
;o  y( 

k 

a  +6 


which,  after  simplification,  results: 


«v  ■  6v  ■  Cos“k  ■  y\  sinak)« 

/j  /0  J  J 


=  "  yj  +  yo  +  xj  Sinotk  +  yj  Cosak 


(B29) 


ARTICULATING  CONSTRAINT  SURFACES 
From  equations  (85)  and  (86) 


yc  =  A1  *  Vc  *  Vc2  *  Vc!  *  ASxc4 


(B30) 


in  equation  (B30)  let 


x  =  x  +  a 
C  C  X. 


Ir 

y  =  y  +6 
7  c  7c  y 


then 
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yc  +  6y  =  A1  +  A2(xc  +  6x  5  +  A3(xc  +  6x  3  + 


*  Vxc  *  4xc)4 


or 


yk  +  6  =  A1  +  A2Xc  +  A26x  +  A3(xc)2+  2 A3x 


*  3A  (x^)2«  ♦  A  lx*)*  *  4A5(xJ)3S 

c 


or 


»yc  -  lA2  *  2Vc  *  3VX£>2  *  4Vxc>Xc 


♦  A2xc  *  A3txc)2  *  A4(xc)3  *  A5(xc)<I  '  y 


Similarly,  from  equations  (85)  and  (87) 


y*  =  A,'  +  Alx'  +  Aix'^ 
c  1  2  c  3  c 


in  equation  (B32)  let: 


.  =  x»k  + 


xc  =  xc 


«x» 

xc 


y»  s  y 

7  C  7  C 


I  k  +  * 

V 


which,  after  simplification,  results 


6y'  -  (A;  +  2Aix'k)6  ,  =  A'  +  A'x'k+  A'x'k' 


O 


k 

1 

»  ■ 

i  ' 

»  , 

1:, 

I 

I 

J 


APPENDIX  C 

COMPUTER  PROGRAM,  JNTMDL 

The  following  pages  contain  a  listing  of  the  computer 
program,  JNTMDL,  which  follows  the  numerical  procedures  dis¬ 
cussed  in  this  report.  JNTMDL  was  developed  for  the  two-dimensional 
dynamic  model  of  the  knee  joint  and  produced  the  results 
presented  in  its  discussion. 


c 

c  . 

c 

PROGRAM  JNTMDL 

C 

c  purpose: 

C  Dynamic  analysis  of  a  two-body-sedaented  aodel  of  the  knee 
C  Joint*  in  two  dimensions*  in  response  to  various  dynamic 

C  loadina  functions. 

C 

c  usage: 

C  Porcine  function  amplitude  and  pulse  duration  aust  be  specified. 
C  Time  increment*  delta  T*  may  also  be  varied. 

C  Initial  contact  points  between  the  tibia  and  femur  must 

C  be  specified. 

C 

c  description: 

C  See  DYNAMIC  SIMULATION  OF  THE  ARTICULATING  JOINTS 

C 

c  remarks: 

C  None. 

C 

C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED! 

C  SIHQ  -  solve  a  set  of  simultaneous  linear  eeuations 
C  Various  standard  FORTRAN  library  functions. 

C 

c  method: 

C  Newton-Rsphson  iteration  usind  Newaark  differential 

C  approximations. 


-153- 


ooo  non  nnn  non  oooooorio 


C 

author: 

Mann sour  H*  Moeinzadeh 
DATE? 

October  1980 


IMPLICIT  REAL*8(A-Z) 

INTEGER  IrJ'K'NUM'ITMAXiIT'KS 
DIMENSION  AK(22»22)»D(22)»DELTA(22> 

KS=0 

TIME  INCREMENT  (SEC) 

T=0.0D0 
DELT=0«0001D0 

PULSE  DURATION  (SEC)  AND  AMPLITUDE  (N) 

TPULSE=0.050D0 
AMP=20.0D0 
FI~3»1415?27D0 
FEXTY=0.0D0 
FEXTX=O.ODO 
SAVE=0.0D0 
ALFSAV=300*0D0 

MAXIMUM  SPECIFIED  ITERATION  NUMBER 
ITMAX=100 

COORDINATES  OF  LIGAMENT  INSERTIONS  (M) 

XPl=0«8D-2 
YPl=16.3D-2 
XP2=2.5D-2 
YP2=17.8D-2 
XP3=2.50D-2 
YP3=20,8D-2 
XP4=-0.5D-2 
YP4=21 .30D-2 
C 

C  COORDINATES  OF  LIGAMENT  ORIGINS  (M) 


N 1 X=-GAMA*P/DSQRT < 1 +P*P ) 

N 1 Y =G AM A/DSQR  T ( 1 +P  *P ) 

N2XP=-BETA*Q/DSQRT< 1+Q*Q> 

N2YP=8£TA/BSQRT (1+0*0) 

INITIAL  FLEXION  ANGLE 

ALFA=BATAN<<P-Q)/(1+P*Q) ) 

ALFA-ALFA+F‘1 

SN=HSIN(ALFA) 

CN:BCCS(ALFA) 

TIBIAL  CENTER  OF  MASS 

XO=XC-XPC*CN+YPC*SN 
YO- YC -XPCtSN- YPCtCN 
FN-O.ODO 

CQ£FCICIENT  OF  FRICTION 
MEli-O.ODO 

INITIAL  LINEAR  AND  ANGULAR  VELOCITY  AND  ACCELERATION  OF 
TIBIAL  CENTER  OF  MASS 

XOMl=XO 

XODMl=O.ODO 

XODDMl=O.ODO 

YOMl=YO 

YODH1=0»ODO 

YODDM1=O.ODO 

ALFM1=ALFA 

ALFDM1-O.ODO 

ALFDDl=O.ODO 

X1=X0+XP1*CN-YP1*SN 
X2=X0+XP2*CN-YP2*SN 
X  3‘ X0+XP3*CN- YP3HSN 
X4=X0+XP4*CN-YP4*SN 

Yl=YO+XPUSN+YPl*CN 

Y2=Y0+XP2*SN+YP2*CN 

Y3=Y0+XP3*SN+YP3*CN 

Y4=Y0+XP4*SN+YP4*CN 


non  o  no 


C 


101  CONTINUE 

PF’RIM=2*C1+6*B1*XC+12*E1*XC**2 

QF'R1M=2*C2 

GAMA-PF'RIM/DABS  ( PPRIN ) 
BETA=GPRIM/DABS(QPRIH) 

C 

INITIAL  LENGTHS  OF  LIGAMENTS 

L1=DSQRT  ( (X5-X1 )**2+< Y5-Y1 >tt2> 
L2=DSQRT  <  ( X6-X2 )  «2+  ( Y6-Y2 )  «2 ) 
L3=DSGRT  <  ( X7-X3 )  W2+ <  Y7- Y3 )  M2 ) 
L4=BSQRT ( ( X8-X4 ) M2+ ( Y8-Y4 >M2 ) 

IF ( T .EQ.0.0)L1I=L1 
IF ( T .EQ.O.O)L2I=L2 
If  :T.EQ.0.0)L3I=L3 
IF(T.EQ.0.0)L4I=L4 

LIGAMENT  FORCES 

IF  (LI  *GT  .L1I)  ABSF1=K1*(L1-L1I)M2 
IF(L2.GT.L2I)  ABSF2=K2* ( L2-L2I > M2 
IF(L3.GT.L3I)  ABSF3-K3* ( L3-L3I > M2 
IF(L4.GT.L4I)  ABSF4*K4*(L4-L4I)M2 
C 

IF <L1«LE«L1I)  ABSF1*-ABSF1 
IF<L2.LE.L2I>  ABSF2=-ABSF2 
1F(L3.LE.L3I)  ABSF3=-ABSF3 
IF(L4.LE.L4I)  ABSF4*-ABSF4 


F1X=ABSF1*(X5-X1)/L1 
F2X= ABSF  24  <  X&-X2 ) /L2 
F3X~ ABSF34 ( X7-X3 ) /L3 
F4X=ABSF4*<X8-X4)/L4 
F1Y- ABSFltl Y5-Y1 )/Ll 
F2Y=ABSF2*<Y6-Y2)/L2 
F  3Y- ABSF34 ( Y7-Y3) /L3 
F4Y-A8SF44(Y8-Y4)/L4 

ALFDEG=ALFA*180 . 0/PI 
IF<T *LT.0.05)G0  TO  8 
IF(ALFSAV  .LT.  ALFDEG)  SOTO  999 
8  CONTINUE 
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ALFSAV=ALFDEG 


C 

t  MOMENT  ARM 

C 

M0MARM=0<0D0 

C 

C  TOTAL  EXTERNAL  ENERGY  <N-H) 

C 

EXTENR=DABS<  XO-TEHPXO ) tFEXTX+DABS ( YO-TEMPYO ) *FEXTY+ 
♦ ( FEXTXtMOMARM ) *DABS < TEMALF-ALFA > 

C 

TOTEXT=EXTENR+SAVE 

C 

SAVE=TOTEXT 

TEMPXO=XO 

TEMPYO=YO 

TEMALF=ALFA 

C 

C  TOTAL  INTERNAL  ENERGY 

C 

P0TENl=<Kl*(Ll-LlI)M3>/3 
POTEN2=  ( K2*  ( L2-L2I )  M3)  /3 
P0TEN3=(K3*(L3-L3I)**3)/3 
F'0TEN4= ( K4*  ( L4-L4I  >**3  >/3 
C 

IF<Ll.LE.LlI)POTENl=O.ODO 
IF(L2.LE.L2I) P0TEN2=O . ODO 
IF<L3<LE<L3I ) POTEN3=0 . ODO 
IF(L4.LE.L4I)P0TEN4=0.0D0 


POTENR=POTEN 1 +POTEN2+POTEN3+POTEN4 

KINENR= <  MOMINRIALFDM1M2 ) /2+SEGMAS*  <  X0DM1«2+Y0DM1**2 ) /2 

TOTINT=POTENR+KINENR 

YC=Aim*XC+Cl*XCU2+Dl*XC**3+El*XC*«4 

YPC=A2+B2*XPC+C2*XPC**2 

IF(T.GT.O.O)  GO  TO  102 
URITE<6*15) 

15  FORMAT < '  '»3X»'T'»4Xf ,ALFA'f5X»,ALFA<'»T28»'X0'»T38»,Y0'»9Xf 'X»' 
♦*9X*,Y<'*T68*'XC'*T78*'XPC'*T88*'FN'*T95*/MC<LIG'*T103*'LC<LIG'* 
♦Till* 'AC<LI6' *T119* 'PC<LIG'*T128» 'IER<#' *2X* 'EXT-ENR'*2X* 'INT-ENR 
♦' *6X* ' YC' *7X* 'YPC' *//) 

102  WRITE(6«16)T* ALFDEG*ALFDD1»X0*Y0*X0DDH1*Y0DDM1*XC*XPC*FN»ABSF1* 


non 


IABSF2 , ABSF3  > ABSF4 » IT* TOTEXT , TQTINT » YC  >  YPC 
16  FORHAT ( '  ' *F7.5*T10»F6.2i T16fF8.1»T26>F7*4»T36»F7.4>T45»F8.2fT56» 
iF8.2*T66>F8*5>T76*F7.5>T85»F8.2»T94>F7.2»T102*F7.2>T110»F7.2t T118 
$»F7.2»T128» I4»3X»F7.4»2X»F7.4*3X»F8*5»3X»F8.5) 

C 

TZERO=TPULSE 

APPLIED  DYNAMIC  FUNCTION 

PULSE =AMP*EXP ( -4 . 73*  ( T/TZERO) **2  >*DSIN (PI*T/TZER0 ) 

L 

FEXTX=PULSE*DCOS(ALFA-PI ) 

FEXT Y=PULSE*DSIN ( ALFA-PI ) 

C 

IF  <  T . GE . BELT .  AND . T . LE . TPULSE )FEXTX=FEXTX 
IF  <  T . GE . BELT . AND . T  *LE . TPULSE ) FEXTY=FEXTY 
I F  •  T . GT . TPULSE ) FEXTX=0 . ODO 
IF ( T . GT . TPULSE ) FEXT Y=0 . ODO 
C 

T=T+DELT 

IT=0.0 

C 

C  NEWION-RAPHSON  ITERATION  PROCESS 

L- 

1  CONTINUE 
IT::IT+1 

IF< IT.GT . ITMAX)WRITE<6f 39) 

39  FORMA '  >T2r  'NUMBER  OF  ITERATIONS  EXCEEDED  THE  SPECIFIED 

t  ITERATION  NUMBER  !  EXECUTION  ABORTED*.* . ') 

IF(IT.GT.ITMAX)GO  TO  999 
C 

SN=DSIN(ALFA) 

CN=DCOS(ALFA) 

C 

C  INITIALIZATION  OF  CAK3  AND  CD]  MATRICES 
C 

DO  5  1=1 t 22 
D(I )=0.0 
DO  5  J=1 »22 
AK< I > J)=0*0 
3  CONTINUE 

r* 

L 

C  COMPONENTS  OF  CD3  MATRIX 
C 

n < 1 > =-Xl+XO+XPl*CN-YPl*SN 
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D < 2  > =-X2+X0+XP2*CN-YF'2*SN 
D  <  3 ) =-X3+X0f  XP3ICN- YP31SN 
D ( 4 ) =-X4+X0+XP4*CN-YP4*SN 

w 

Ii..5)=-Yl+Y0+XF1*SN+YP1*CN 
£i  <  6 )  =-Y2+Y0+XP2KSN+YP2tCN 
D(7)  »~Y3+Y  0PXF'3*SN+YP3*CN 
D ( 8 ) --Y4PY0tXP4*3N+YP4<CN 
C 

[i  (9)=-  SE6MAS*  ( ( 4 . 0/DELTW2 ) * ( X0-X0M1 ) - <  4 . 0/DELT )  *XODHl-XODDMl ) 

*- ( F1X+F2X+F3XPF4X ) +FN* ( MEU*N1 Y-N1X  > *GAHA-FEXTX 

C 

[i  (10 ) -SCuHASt ( <  4 . 0/DELTU2  > * ( Y0-Y0H1 ) - ( 4 . 0/DELT ) HY0DM1-Y0DDM1 ) 

$ - ( FI Y+F2Y+F3Y+F4Y ) -FN* ( MEU*N1X+N1 Y  >*GAMA-FEXTY 

C 

Dai)=-((Xl-XO)*FlY+(X2-XO)*F2Y+(X3-XO)*F3Y+(X4-XO)*F4Y-(Yl-YO)* 
*F1X-(Y2-YO>*F2X-<Y3-YO>*F3X-<Y4-YO>*F4X>-FN*GAMA*((XC-XO)*<N1Y+ 
$MEU*N1X)-<YC-Y0)*<N1X-MEU*N1Y>)+MQMINR*(<4.0/DELT**2)*<ALFA-ALFM1 
$ ) - ( 4 . 0/DELT ) *ALFDM1-ALFDB1 ) +FEXTX*MOMARM 
D(12)=X0-XC+XPC*CN-YPC*SN 
D(13)=Y0-YC+XPC*SN+YPC*CN 
IK15)-Al+Bl*XC+Cl*XCi|c*2+DUXC«3+El*XC**4-YC 
D  ( 16 )  =A2+B2*XF'C+C2*XPC**2-YPC 
Ii  ( 1 7 ) =B1 +2*C1 *XCf 3*U1  *XC**2t4*El *XC**3-P 
D  a  8 )  ( P/DSQRT  ( 1 +P**2 ) )  *GAMA-N1 X 

D(  19 )  =  ( 1 . 0/DSQRT  <  1+P«2 ) )  tGAMA-NlY 
K20)=B2t2*C2*XPC-Q 
H ( 21 ) = ( -Q/DSQRT < 1 +0**2 ) ) *BETA-N2XP 
3 i 22 ) = (1 . O/DSQRT ( 1 +Q**2 ) >*BETA-N2YP 
C 

C  COMPONENTS  CiF  CAK3  MATRIX 
C 

AK< 1 ; 1 ) --1 . ODO 
AKa»3)=XPl*SNtYPl*CN 
AK  1 » ? ) - 1 . ODO 
t 

AKv2>l)--1.0D0 
Ai<(2»3>=XP2*3N+YF'2*CN 
AK(2>8)~1 .CDO 
C 

AK(3»1)=-1.0D0 

AK(3;3)-XP3*SNtYP3tCN 


AK(4»i)“-l«0D0 
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fth : } j  3  > -XP4#SN+YP4*CN 
nK  •:  ‘It  10  >  =  1  »0D0 
C 

hK(5»2)==-1.0D0 

hK(A»2) -~1*0D0 
AK<7»2)=-1.0DO 
AK(8>2)=-1»0D0 

r' 

AK ( 5 » 3 ) ( XP1 tCN-YPl<SN ) 

AK(6>3)--( XP2*CN- YP2*SN ) 

AK ( 7 i 3 )-- <XP3*CN-YP3tSN) 

AK ( 8  >  3 ) *- < XP4*CN- YP4*SN) 

C 

AK(5»11)=1.0D0 
AK(A»12)=1.0D0 
AK(7»13)=1*0D0 
AK(8> 14)=1 »0D0 
C 

AK(9  ,4>=-<MEU*NlY-NlX)*GAHA 
AK(9  rl6)=-MEU*FN*GAMA 
AK(9  >15)-FN$GAMA 
AK(9  r  1  >  =-4 .0ISEGHAS/DELTM2 

C 

AK(10»4)=(NEU*N1X+N1Y)*GAMA 
AK ( 1 0 » 1 6 ) =FN*GAMA 
AK'10,15)=MEU*FN*GAMA 
AK  ( 1 0  >  2 )  =-4 . 0*SEGHAS/DELT«2 
C 

AK(llf7)=FlY 

AK(11»8)=F2Y 

AK(11»?)=F3Y 

AK(ll»10)=F4Y 

AK  ( 1 1 . 1 )  =-  ( ( F 1 Y+F2YFF3  Y+F4Y )  +FN*  ( HI  Y+MEMN1 X )  *GAMA ) 

AK(11»11)=-F1X 

AK(llrl2)=-F2X 

AK<11»13)=-F3X 

AK(11»14)=-F4X 

AK (1 1  *  2 ) = < ( F 1 X4F2X+F3X+F4X) +FN*  <  MIX-HEUtNl Y ) SGAHA ) 
AK( 11 »4)=(XC-X0)t<NlYFMEU$NlX)-(YC-Y0)t<NlX-MEU<NlY) 
AK<11»15)=( HEUtFNt ( XC-XO ) -FH* ( YC- YO ) ) *GAMA 
AK (1 1 » 1 6 > = ( FN* ( XC-XO  >  FHEUIFN* < YC-YO ) ) IGAHA 
AK(ll»5)*FHt(HlY+HEU<tNlX)$GAMA 
AK  <  1 1  >  1 9  > =-FN* ( N 1 X -NEU*N 1 Y ) *G AH A 
AK  <  1 1 » 3 )  =-4 . 0IH0HINR/DELTM2 
C 
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AK(12»5)=1.0D0 
AK(12»1)*-1.0D0 
AK ( 12  r  3 ) *XPC4SN+ YPCtCN 
AK(12»6)=~CN 
AK(12-20)=SN 
C 

AK(13t2)=-1.0D0 
AK ( 13 • 3) = YPCtSN-XPCtCN 
AK(13rl9)-1.0D0 
AKC13»6>=-SN 
AK(13f20)>-CN 
C 

AK ( 14 1 22)*P*SN+CN 
AK(14,21)=Q*SN-CN 
AK ( 1 A t 3) =CN* ( 1+PtQ ) +SN* ( P-Q) 

C 

AK( 15» 19)=1 .ODO 
AK(15»5)*-P 
C 

AK(16»20)=l*0DO 

AK(16»6)=-Q 

C 

AK(17»21)=1.0D0 

AK (1 7  t  5 ) =- < 2*C1+4*D1*XC+12«E1*XCM2 ) 

C 

AK(18»15)=1.0D0 

AK(18*21)*GAMA/<1+P*P>«1.5 

C 

AK(19f 1A)=1*0D0 
AK( 19*21 )=BAMAtP/(l+PtP)ttl.5 
C 

AK(20»22)al*0D0 

AK<20*A)=-2*C2 

C 

AK(21»17)*I.0D0 

AK(21»22>CBETA/(1+QM)W1*5 

C 

AK(22rl8)=1.0DO 
AK ( 22 * 22 ) »8ETA*Q/ ( 1 +8*0 > *«1 . 5 

CALCULATION  OF  COMPONENTS  OF  DELTA  MATRIX 

CALL  SIMQ<AK*D*22*KS) 

C 

DO  44  1*1*22 
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0ELTA<I)=B<I> 
44  CONTINUE 


IF(KS.EQ.l)  URITE<6»3) 

3  FORMAT ( '  '*'KS=1  , SINGULAR  MATRIX  • EXECUTION  ABORTED  . ') 

IF(KS.EQ.l)  GO  TO  999 

NEW  VALUES 


X0=X0+DELTA<1> 

Y0=Y0+DELTA<2) 

ALF  A=ALFA+DELT  A ( 3 ) 

FN=FN+D£LTA<4) 

XC=XC+D£LTA(5) 

XPC*XPC+D£LTA ( 6 ) 

X1=X1+0ELTA(7) 

X2=X2+DELTA<8> 

X3=X3+DELTA(9) 

X4=X4+DELTA(10) 

Y1-Y1+DELTA(11) 

Y2--Y2FDELTAU2) 

Y3; Y3+DELTA(13) 

Y4=Y4+DELTA( 14) 

N1X=N1X+DELTA(15) 

N1Y-N1YFDELTA<16) 

N2XF-N2XP+DELTA( 17) 

M2YP=N2YP+DELTA( 18) 

YC=YCfDELTA(19) 

YPC=YPCFD£LTA<20) 

P=P+DELTA(21) 

Q=Q+DELTA(22) 

CONVERGENCE  TESTS 

0MEG=0.0001D0 


DO  555  1=1 »22 

IF ( DABS<  DELTA( I ) ) «GT » OHEG )  GO  TO  1 
555  CONTINUE 

NEUMARK  DIFFERENTIAL  APPROXIMATIONS 

CALCULATION  OF  VELOCITIES  AND  ACCELERATIONS  AT  TIME  T 

X0DDT= ( 4 . O/DEL Tt *2 ) * ( XO-XOM1 ) - ( 4 • O/DELT ) *X0DM 1 -XODDM 1 


XODT =XODMl + ( DELT/2 . 0 ) tXODDMlt ( DELT/2 . 0 ) tXODDT 
C 

XOHl=XO 

XODHl=XODT 

XODDMl=XODDT 

C 

YOHDT=  ( A .  0/DELTM2  >  *  ( YO-YOM1 ) -  ( 4 . O/DELT )  *YODMl- YODDM1 
YODT= Y0DH1 +  <  DELT/2 . 0  > *YODDMl + <  DELT/2 . 0 ) *YODBT 
C 

YOMl=YO 

YODMl=YODT 

YOHDhl=YODDT 

C 

ALFDDT=  ( A .  0/DELTM2)  *  ( ALF  A-ALFM1 )  -  ( 4 . 0/DELT )  *ALFDM1-ALFDD1 
ALFDT= ALFDM 1  +  <  DELT/2 . 0 ) tALFDD 1 + ( DELT/2 . 0 ) 4ALFDDT 

C 

ALFM1=ALFA 

ALF0M1=ALFBT 

ALFDD1=ALFDDT 

C 

GO  TO  101 

c 

999  STOP 
EMD 
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SUBROUTINE  SIMQ(A*B*N»KS> 

C 

C  . . 

c 

purpose: 

OBTAIN  SOLUTION  OF  A  SET  OF  SIMULTANEOUS  LINEAR  EQUATIONS* 
AX=B 


usage: 

CALL  SIMQ(A*B*N»KS) 

DESCRIPTION  OF  PARAMETERS: 

A  -  MATRIX  OF  COEFFICIENTS  STORED  COLUMNWISE*  THESE  ARE 
DESTROYED  IN  THE  COMPUTATION.  THE  SIZE  OF  MATRIX  A  IS 
N  BY  N* 

B  -  VECTOR  OF  ORIGINAL  CONSTANTS  (LENGTH  N).  THESE  ARE 
REPLACED  BY  FINAL  SOLUTION  VALUES*  VECTOR  X. 

N  -  NUMBER  OF  EQUATIONS  AND  VARIABLES.  N  MUST  BE  .GT.  ONE. 
KS  -  OUTPUT  DIGIT 

0  FOR  A  NORMAL  SOLUTION 
1  FOR  A  SINGULAR  SET  OF  EQUATIONS 


remarks: 

HATRIX  A  MUST  BE  GENERAL. 

IF  MATRIX  IS  SINGULAR*  SOLUTION  VALUES  ARE  MEANINGLESS. 

AN  ALTERNATIVE  SOLUTION  MAY  BE  OBTAINED  BY  USING  MATRIX 
INVERSION  (MINV)  AND  MATRIX  PRODUCT  (GMPRD). 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED: 

NONE 

method: 

METHOD  OF  SOLUTION  IS  BY  ELIMINATION  USING  LARGEST  PIVOTAL 
DIVISOR.  EACH  STAGE  OF  ELIMINATION  CONSISTS  OF  INTERCHANGING 
ROWS  WHEN  NECESSARY  TO  AVOID  DIVISION  BY  ZERO  OR  SMALL 
ELEMENTS. 

THE  FORWARD  SOLUTION  TO  OBTAIN  VARIABLE  N  IS  DONE  IN 
N  STAGES.  THE  BACK  SOLUTION  FOR  THE  OTHER  VARIABLES  IS 
CALCULATED  BY  SUCCESSIVE  SUBSTITUTIONS.  FINAL  SOLUTION 
VALUES  ARE  DEVELOPED  IN  VECTOR  B*  WITH  VARIABLE  1  IN  B(l>* 

VARIABLE  2  IN  B(2>*... . *  VARIABLE  N  IN  B<N). 

IF  NO  PIVOT  CAN  BE  FOUND  EXCEEDING  A  TOLERANCE  OF  O.O* 

THE  MATRIX  IS  CONSIDERED  SINGULAR  AND  KS  IS  SET  TO  1.  THIS 
C  TOLERANCE  CAN  BE  MODIFIED  BY  REPLACING  THE  FIRST  STATEMENT. 
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IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  A<l)rB(l) 

FORWARD  SOLUTION 

T0L=0*0 

KS=0 

JJ=-N 

DO  A5  J=1»N 

JY=J+1 

JJ=JJ+N+1 

BIGA-0 

IT=JJ-J 

DO  30  I-J>N 


SEARCH  FOR  MAXIMUM  COEFFICIENT  IN  COLUMN 
1J=IT+I 

IF(DABS(BI6A)-DABS(A(IJ) ) )  20>30f30 
20  BIGA=A(IJ) 

IMAX=I 
30  CONTINUE 

TEST  FOR  PIVOT  LESS  THAN  TOLERANCE  (SINGULAR  MATRIX) 

IF(DABS(BIGA)-TOL)  35»35»40 
35  KS=1 
RETURN 

INTERCHANGE  ROWS  IF  NECESSARY 

40  Il*J4N*(J-2) 

IT’IMAX-J- 
DO  50  K*JiN 
1 1*1 1-fN 
I2=I1+IT 
SAVE»A(I1) 

AC  II )=A< I2> 

A(I2)«SAVE 

DIVIDE  EQUATION  IY  LEADING  COEFFICIENT 
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SO  A<Il)=A(ll)/BIGA 
SAVE=B( I MAX) 

B(IMAX)=B<J) 

B( J)=SAVE/BIGA 

ELIMINATE  NEXT  VARIABLE 

IF(J-N)  55f70»55 
55  IQS=N*<J-1> 

DO  65  IXeJYfN 

IXJ=IQS+IX 

IT=J-IX 

DO  60  JX=JY»N 

JXJX=N*(JX-1)«X 

JJX-IXJX+IT 

60  A( IXJX)=A< IXJX)-(A(IXJ)tA( JJX) ) 
65  B(IX)sB(IX)-(B<J)tA(IXJ)> 

BACK  SOLUTION 

70  NY=N-1 
IT=N«N 

DO  80  JS1*NY 
IA=IT-J 
1B*N-J 
IC=N 

HO  80  K=l» J 

B(IB)=B(IB)-A(IA)*B(IC) 

IA=IA-N 
30  IC=IC-1 
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